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ABSTRACT 

A second-order epsilon method Is developed for trajec- 
tory optimization problems. The method is applied to sever- 
al aircraft and missile performance and air combat maneuver- 
ing problems. Heavy emphasis is placed on the realistic 
modeling of the flight vehicle's motion and maneuvering 
limitations . 

The proposed optimization technique, which is an exten- 
sion of Balakrishnan’ s epsilon method, uses either the full 

second-order Newton-Raphson method or the "modified" Newton- 

/ 

Raphson method to minimize the epsilon functional. The full 
Newton-Raphson method exhibits terminal convergence charac- 
teristics superior to the "modified" method, whereas the 
"modified" method is generally superior in the initial 
stages of a problem. An algorithm is developed which uses 
both techniques in a complementary way. 

A new penalty functional which has desirable theoretical 
properties and exhibits excellent computational behavior is 
introduced to treat state and control inequality constraints. 
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I . INTRODUCTION 



The objective of the research reported on herein Is to 
develop a method of solving realistic problems In aircraft 
and missile performance optimization. Optimization problems 
of this type have been the subject of considerable research 
[Refs. 1 , 2 , 3j 5, 6, and ?]• The mathematical models 
used In these references are the products of many simplifi- 
cations and assumptions. Typically, the degree of simplifi- 
cation used to render these problems solvable by some opti- 
mization technique Is such that the solutions obtained are 
of limited practical value. This Is particularly true In 
the modeling of aircraft maneuvering limitations, such as 
aerodynamic stall, maximum structural load factor, and 
placard Mach number, which require the use of multiple 
state and control Inequality constraints. Since these 
limitations play an extremely Important role In maneuvering 
flight and air combat, they must be modeled as accurately 
as possible. It Is, therefore. Imperative that the optimi- 
zation technique used to solve the problems posed herein 
be capable of handling state and control inequality 
constraints with relative ease. 

Balakrlshnan’ s epsilon method is an attractive optimiza- 
tion technique because of the natural manner In which state 
and control Inequality constraints are introduced. The 
epsilon method Is a penalty method In which terms are added 
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to the performance measure to penalize deviations from the 
state equations written as equality constraints. Likewise, 
state and control inequality constraints may be treated by 
the addition of appropriate penalty terms to the performance 
measure. The resulting augmented performance measure is 
minimized by an appropriate algorithm for solving unconstrained 
optimization problems. 

The optimization technique used most successfully in 
the literature with the epsilon method is a "modified" 
Newton-Raphson technique, hereafter referred to as the MNR 
technique. In this method certain second-order terms pre- 
sent in the full Newton-Raphson formulation are neglected. 

It is argued [Ref. 8] that since computer storage and time 
requirements to compute these terms are large, and since 
satisfactory results can be obtained over a large class of 
problems without the terms, their inclusion is not Justified. 
For these reasons the full Newton-Raphson formulation, here- 
after referred to as the FNR technique, has not been 
previously utilized with the epsilon method. 

Difficulties were experienced by the author, however, 
in applying the MNR technique to problems of the type formu- 
lated in this dissertation. The MNR technique was not 
effective in problems with state equations and multiple 
inequality constraints resulting from a realistic modeling 
of the flight vehicle's motion and maneuvering limitations. 

For this reason the FNR method was investigated and 
found to be feasible in terms of computational storage and 
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time requirements. The FNR method exhibits terminal con- 
vergence characteristics superior to the MNR method although 
the MNR method is generally superior in starting a problem. 
Problems not solvable by the MNR method alone were solved 
by an algorithm which uses both methods in a complementary 
way. The power of the FNR technique close to the minimum 
can also be used to advantage to obtain a family of optimal 
trajectories for different end conditions. The optimal 
trajectory for one set of end conditions is used as a first 
guess for the optimal trajectory for a neighboring set of 
end conditions. 

Several simplified problems in aircraft performance 
optimization were attempted initially to gain experience 
with the epsilon method. In these problems inequality 
constraints were treated by using interior penalty func- 
tionals of the type recommended by Fiacco, McCormick, and 
Jones [Refs. 9 and 10]. Computational results were unsatis- 
factory. Difficulties were experienced in keeping the con- 
strained state or control completely admissible; a require- 
ment for the success of an algorithm with this type of 
penalty functional. To alleviate this difficulty a new 
penalty functional for inequality constraints is introduced 
which exhibits excellent computational behavior. The pro- 
posed functional has performed well in computation with up 
to eight inequality constraints represented in a single 
problem. 
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The thesis is divided into eight sections. In Section 
II the epsilon method is presented. The FNR and MNR tech- 
niques are derived and discussed. The effectiveness of 
the FNR method as opposed to the MNR method is demonstrated 
by a scalar example. Finally, the computational experience 
gained with both methods in solving realistic performance 
problems is presented. In Section III the method of treating 
state and control Inequality constraints is presented. The 
author’s experience with interior penalty methods is related 
and a new penalty functional is proposed. Computational 
experience with the new penalty functional is related and, 
finally, several desirable theoretical properties of the 
new penalty functional are presented. In Section IV the 
algorithm developed for minimizing the epsilon functional 
by either the MNR or FNR methods is presented. In Sections 
V, VI, and VII three aircraft and missile performance opti- 
mization problems are solved. These problems are pertinent 
and realistic in their operational applicability. The 
three-degree-of- freedom models are the same as those used 
in basic aircraft performance analysis. Finally, the 
summary and conclusions are presented in Section VIII. 



\ 



9 



II. THE EPSILON METHOD 



This section describes the epsilon method and reviews 
the significant contributions of other Investigators. The 
full Newton-Raphson (FNR) equations for minimizing the aug- 
mented performance measure are derived and compared to the 
modified Newton-Raphson (MNR) equations published elsewhere 
[Refs. 8, 11, and 12]. 

A. DESCRIPTION OF THE EPSILON METHOD 
1 . Statement of the Problem 

A dynamic system characterized by the nonlinear 
state equations 



x(t) = f[x(t) ,u(t) ,t] (2.1) 

Is to be controlled to minimize the performance measure 

J(x,u) = h[x(T),T] + / g[x(t ) ,u(t) , t]dt (2.2) 

- ~ 

where x(t) Is an n x 1 state vector and u(t) Is an Jl x 1 
control vector. State and control Inequality constraints 
are omitted for the present. In Section III the Inclusion 
of these constraints Is discussed In detail. 

2 . The Augmented Performance Measure 

In the epsilon method as proposed by Balakrlshnan 
[Refs. 13 and 1^], the performance measure (2.2) Is augmented 
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by a penalty functional which involves a weighted integral 
of the Euclidean norm of the state equations written as 
equality constraints. The augmented performance measure is 




o 



= J(x,u) + ^ J (x,u) . 



(2.H) 



The weighting factor e is a positive quantity. 



3 . Behavior as e -*• 0 



As e is reduced, the penalty term J is more heavily 

s 



weighted, thereby placing greater emphasis on satisfying 
the state equations. Balakrishnan [Refs. 13 and 1^] and 
Taylor [Ref. 11] have shown that under appropriate assump- 
tions as e -»■ 0, the epsilon method yields the necessary 
conditions of optimality obtained by applying Pontryagin's 
minimum principle [Refs. 15 and l6]: 



x*(t) = [x*(t),u*(t),t] 



(2.5) 




( 2 . 6 ) 



and 



H[x^(t),u*(t),p*(t) ,t] <. H[x*(t),u(t),p*(t),t] (2.7) 



where H is the Hamiltonian function defined as 



\ 
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H[x(t ) ,u(t) ,p(t) ,t ] = g[x(t) ,u(t) ,t ] + p'^(t)f[x(t) ,u(t) ,t], 

( 2 . 8 ) 

p(t) is the costate or adjoint vector, u*(t) is an extremal 
control vector, and x*(t) is an extremal trajectory. The 
assumptions made are that the minimization problem has a 
unique solution with x(t) absolutely continuous for each e, 
and that f and g are continuously differentiable in x and 
u [Ref. 11]. Thus, under appropriate assumptions, it can 
be shown that as e -»■ 0, the epsilon method yields the 
results of Pontryagin's minimum principle. That is, if the 
optimal control u*(t,e) of equation (2.3) exists for each 
e, that solution will approach the optimal control u*(t) 
of equation (2.2) as e 0. 

. State Equation Integration 

It should be noted that the epsilon method is a 
non-dynamic method in that the state equations are not 
integrated during the minimization process. Once the aug- 
mented performance measure has been minimized a check on 
the degree of satisfaction of the state equations can be 
obtained by integrating the state equations with the optimal 
control . 

V 

B. MINIMIZING THE AUGMENTED PERFORMANCE MEASURE 
1 . Sequence of Unconstrained Problems 

Once the augmented performance measure (2.3) is 
formulated, any unconstrained optimization algorithm can 
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be applied to it. A sequence of unconstrained problems 
referred to as sub-problems is solved. In each sub-problem 
e is held constant and a minimization is performed until 
some stopping criterion is satisfied. At this point e is 
decreased and a new sub-problem is commenced using the 
optimum trajectory found in the previous sub-problem as a 
first guess. In this manner a sequence of sub-problems is 
solved until, if convergence occurs, some overall stopping 
criterion is satisfied. 

2 . Unknowns and Time Discretization 

The states and controls can be approximated by any 
orthogonal expansions. The coefficients in these expansions, 
along with all free end conditions, become the parameters 
or unknowns in the optimization. A functional expansion of 
the form (2.9) is convenient because it is continuous and 
the period can be selected so that the value of the expansion 
is zero at the end points. Since the problems solved involve 
time-invariant systems, t^ is selected as zero and the states 
and controls are written as 



y(t) 



x(t ) 



u(t ) 



y(o) 



+ 

T 



t + D 



sin 



trt 



T 

sin 2trt 
T 



( 2 . 9 ) 



sin Mut 
T 



where M is the number of harmonics used and 
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♦ 









# fi^" 



D 



( 2 . 10 ) 



D 



X 




is an (n + Jl) x M matrix of coefficients. The derivative 
of the expansion of the state vector given in equation (2.9) 
with respect to time is required and is given by 



x(T) - x(0) 

x(t) = =: 

T 



TL 

T 



cos 



irt 

T 



+ D 

..X 




2irt 

T 



( 2 . 11 ) 



Mtt 

T 



cos 



Mirt 

T 



The objective is to find the D matrix along with the values 
of the free end conditions which minimize equation (2.3) 
for a given e. In order to perform this minimization, the 
time interval T is divided into (K - 1) sub-intervals each 
of duration At so that there are K discrete time points. 

The augmented performance measure given by Equation (2.3) 
is written as 



Ja(x,u,e) 




2 

[x^( t)-f^(x,u) ] dt + 




2 

(x,u)] dt 



^T ^ 2 1 (2.12) 

+ ... + [x^(t)-f^(x,ul] dtj 

r'^ 

+ I g[x,u]dt + h[x(T),T]. 

•'O ~ ~ 



li| 



[ 



Suppressing the arguments for clarity, equation (2.12) is 
expanded to yield 



1 r • P /*2At p ^ ivAi/ , 

"^a ‘t[J^ ^V^l) y -"X .... (v^l)^dt 



KAt 
'(K-l)At 



^At p ^2At p K^t 

/ (x„-f,)^dt +/ (x„-f„)^dt + ••• + / (iL-f^) 

^0 ^ ^ -'At ^ ^ *XK-l)At ^ ^ 



'dt 



(2.13) 



^At . ^At „ 

4 (vy <vy it 



+ ••• + 



/. KAt 

/ 

iK-DAt 






/^t /Q.bX> ^ KAt 

/ gAt +/ gAt + ••♦ + / gAt + hCx(T),T] 

*'0 *'At *TK-l)At 



which can be approximated by 



2 2 

Ja “(5^(0) - f3^[x(0),u(0)]} |i^(At) - f3^[x(At),u(At)]} ^ 

+ ••• + |x^([K-l]At) - f^[x([K-l]At),u((K-l)At)]J ^ + 

|x2(0) - f2[x(0),u(0)]j- |x 2 (At) - f 2 [x(At) ,u(At) ]| ^ 

+ .*•-+ |L((K-l)At) - f„[x((K-l)At),u((K-l)At)]) ••• + 

2 '■ ^ Z 

|x^(0) - f^[x(0),u(0)]| |x^(At) - f^[x(At),u(At)]J ^ 

+ ♦♦♦ + |3^((K-l)At) - f^[x((K-l)At),u((K-l)At)]} ^ + 

|g[x(0),u(0)]|At + |g[x(At),u(At)]|At + ...+ 

|g[x( (K-l)At) ,u( (K-l)At) ]jAt + h[x(T) ,T] .' 
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Hence, the augmented performance measure can be written as 



J 



2 



+ w 



+ 



+ w 



2 



(2.15) 



a 



w 



1 



2 



Q 




( 2 . 16 ) 



where w is a Q x 1 column vector. The first K elements of 
w are 



etc. The form of equation (2.16) is convenient for computer 
programming the epsilon method and for the derivation of 
the minimization techniques that follow. In minimum time 
problems where the performance measure is given by 

/•T 

J = / dt. (2.18) 

"b 

one element of w of the form 



is used to represent equation (2.18). In these type problems 
the number of time points K is held constant during the 
minimization in order to keep the dimensions of all vectors 
and matrices constant and the time interval At is minimized. 



w 




Wq = [(K-l)At]^ 



(2.19) 



\ 
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The values of the states and controls required in 
elquatlon (2.1^1) are obtained by evaluating equations (2.9) 
and (2.11) at each time point t = (k - l)At where k = 1,2, 
...,K. Written In discrete form equation (2.9) Is 



y[(K-l)At]-y(0) 

y[(k-l)At] = y(0) + = (k-1) + D 

K-1 



•1n ^(k-1) 

K-1 " 
2Tr(k-l) 



sin 



K-1 



sin 



( 2 . 20 ) 



and equation (2.11) Is 



. x[(K-l)At]-x(0) 

x[(k-l)At] = + D 



(K-l)At 



~x 



TT TT(k-l) 

(K-l)At K-1 

2tt ^ 2TT(k-l) 



(K-l)At K-1 



Mtt MTr(k-l) 

(K-1) At K-1 



( 2 . 21 ) 



- 



A vector of unknowns c Is formed and Is given by 



2yW 







^1,1* *^1,2» •••’ 


*^2,1 


n+Jl,l’ *^n+Jl,2* 





, ^ 2 * * * * ^ At) 



( 2 . 22 ) 
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where ^ is the element In the 1 ^^ row and column 
of p and z^, z^, ...» Zp represent P free end conditions 
some of which occur at t = 0, and others at t = T. Some 
of the Zp's correspond to states and others to controls. 

The last element, At, is present only if time is to be 
minimized. The c vector consists of L elements where 

L = (n+A) X M + P ( 2 . 23 ) 

for all problems except minimum time problems and 

L = (n+Jl) X M + P + 1 (2.24) 

for minimum-time problems. 

With the states and controls given by equation 
( 2 . 20 ) and the augmented performance measure given by 
equation ( 2 . 16 ), the problem has been transformed into a 
parameter optimization problem with the unknowns given by 
c ( 2 . 22 ) . 

3 . Minimization Techniques 

The methods which have received attention in the 
literature for finding c* which minimizes the augmented 
performance measure given in equation (2.3) are the gradient 
method and a "modified” Newton-Raphson method (MNR) . 

The gradient method has been investigated by J. 
Taylor [Refs. 11 and 12] and L. Taylor [Ref. 8] with 
unsatisfactory results. These investigators report that 



in non-linear problems the gradient method frequently 
obtains false minima and requires considerable computation 
time compared to other methods. 

An MNR method in which certain second-order terms 
present in the full Newton-Raphson (PNR) method are 
neglected has enjoyed greater success and requires less 
computation time than the gradient method [Refs. 8, 11, and 
12]. However, in Ref. 8 difficulties are reported with the 
MNR method in non-linear problems. J often begins an 

cL 

oscillation after two or three iterations and does not 
settle to a minimum. In the problems solved herein the 
same oscillations have been observed when the MNR algorithm 
has been used. Convergence to the minimum, when it does 
occur, is typically very slow. Typical performance of the 
MNR method is shown in Figure 1. 




Figure 1 

Augmented performance measure vs. iteration 
number-MNR method 
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The PNR method has not been used by other investi- 
gators with the epsilon method because of the increased 
computer time for each iteration, the additional storage 
space required, and a significantly increased analytic 
workload involved in deriving second partial derivatives. 
Because of the poor performance of the MNR method on problems 
of the type solved herein, the FNR method has been investi- 
gated in detail in this work. With careful programming the 
computer time for each iteration and storage space required 
for the FNR method has been reduced to an extent which makes 
the method computationally feasible. 

^ . The Full and Modified Newton-Raphson Equations 

The FNR equations for finding c* are derived here 
in a manner which permits the MNR equations derived in the 
literature [Refs. 8 and 11] to be obtained by neglecting 
a term. 

The augmented performance measure is expanded in a 
Taylor series Including up to second-order terms and is 
written as 

J (c+Ac) = J (c) .+ (2.25) 

where 
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V J 
c a 



3J. 



3 c, 



and 





3^J 


3^J 


. 3^J 


... 


a 


a 


a 




3Ci^ 


3Ci3c2 






3^J 


3^J . . 


. 9^J 




a 


a 


a 


7 


3 ©2 ^ 


3c2^ 


3c23Cj^ 


c a 










• 


• 


• 




• 


• 


• 




• 


• 


• 




3^J 


3^J . . 


. 3^J 




a 


a 


a 




3Cj^3c^ 


3Cj^3c2 


^ 2 










Solving for the increment 


of J , we have 
a’ 




A 


AJ^ = 


J (c+Ac) - 


J (c) 




a 


a ~ 


a -V 






(V,Ja)/ic 


+ %(Ac)'^(V ) 

O CL 


1 (Ac) . 






( 2 . 26 ) 



( 2 . 27 ) 



( 2 . 28 ) 

( 2 . 29 ) 
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Applying the necessary condition for a minimum, we have 



3(AJ^) 

3(Ac) 



<Wo * Ja>c ^2 = 5 



which when solved for Ac yields 



Ac = -(V^^J^) (V^J^) , 

~ c a c c a c 



If the augmented performance measure is written as 



A T 
= w w 
a ~ ~ 



where w is a Q-dimenslonal column vector, then 



and 



Va “ ’c'H S) 



= 2(V w)^w 

c ^ 



=■ 2V^[(V^w)\] 



= 2 f(V^w)'^(V^w) + (V^^w)\l 

*“ <iv 



The matrix V w is given by 



(2.30) 



(2.31) 



(2.14) 



(2.32) 



(2.33) 



’o!? 



3w^ 


3w2 


3 


3c^ 


3 C 2 


3Cj_ 


3w2 


3w2 


3 Wo 


3 c^ 

• 


3 C 2 

• 

• 


90l 

• 


• 

3“q 


• 

3«Q 


• 

3Wq 


3c^ 


3C2 


9“l 



(2. 3*1) 



V w is a three-dimensional array composed of L matrices each 

~ ~ 3 2wi 

of dimension Q x L which has as its ijk th element ~ 3 ~ c ~ ^ 

dCj^dC^. 

that is 



V 

C 



3^w, 



2 

3 w. 



2 

3 w^. 



3c^^ 


3c^3c2 




3c^3Cj^ 








3^W2 


3^W2 


• • • 


3^W2 


. 






3c^^ 


3Cj^3c2 




3 


3c29Cj^ 






• 

• 


• 

• 




• 

• 


9^W2 


• 




• 


• 




• 


9c29Cl 








9 % 


• • • 


9% 


• 




9^w^ 


3c^^ 


3c^3C2 




3 3 Cj^ 


• 

• 


• 


^ 2 
3 Cl 




9% 


S^w 

^ • • • 


9^Wg 


• 

• 


3^W2 




3C23c^ 


3c2"" 


> 


3c29Cj^ 


• 


9cj_2 






• 


• 


• • • 


• 


• 










o2 

9 Wq 


• • • 










3Cl9Ci 


3Cj^9c2 







(2.35) 
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Substituting equations (2.32) and (2.33) into equation (2.31), 
we have 



Ac =- 






(V w) 

' c~ c 



+ (V. 






V w) 

C-v c 






(2.36) 



This is the full Newton-Raphson equation. The modified 
Newton-Raphson equation can be obtained by neglecting the 
second term in the inverse in equation (2.36), which yields 



Several comments concerning equations (2.36) and (2.37) are 
in order: 

a. the term V w given by equation (2.3^) isaQxL matrix; 

^ m 

b. the term (V w) (V w) is a symmetric L x L matrix; 

C C 

c. the term (V.w) w is an L x 1 vector; 

d. the term V given by equation (2.35) isaQxLxL 

C .V 

three-dimensional array; 

2 T 

e. the term (V^ w) w is a symmetric L x L matrix; 

pm ~ 

f. (V w) is defined as the three-dimensional array 

0 ^ 

obtained by .transposing each individual Q x L matrix 
given in equation (2.35); 

2 T 

g. the result of the operation (V^ w) w is defined as 

an L X L matrix in which the 1^ column is the 
product [1^ w • 



2k 
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5 . A Scalar Illustration of the MNR and FNR Methods 
The potential importance of the second-order term 
neglected in the MNR equations can be illustrated with a 
simple scalar problem in function minimization. The Newton- 
Raphson equation to minimize a function f(x) takes the well 
known form 



which is the form of equation (2.l6) with w taken as a 
scalar, then 




( 2 . 38 ) 



If 



f(x) = w^(x) 



(2.39) 



f * (x) = 2w(x) ^ (x) 



(2.40) 



and 



f"(x) = 2 




^ + 2w(x) 




(2.41) 



The FNR equation (2.36) is 



Ax 




(2.42) 



whereas the MNR equation (2.37) is 



Ax 



(2.i|3) 



-w(x) 



dw 

dx 



(x) 




w(x) 

^ (X) 

dx 



(2.44) 



Applying equations (2.42) and (2.44) to the function 



f.(x) = (x-1)^ + 1 



(2.45) 



In which 



w(x) = [(x-1)^ + 1]^ 



(2.46) 



we obtain 



Ax = - ^ (x-1) 



(2.47) 



for the PNR algorithm and 



Ax 



(x-1)^ + 1 
2(x-l)^ 



(2.48) 



for the MNR algorithm. Tables 1 and 2 show the first few 
Iterations by both methods from an Initial guess of x(0) = 3* 
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MNR Equation ( 2 . 48 ) 


Iteration 


X 


Ax 


f(x) 


1 


3.000 


-1.062 


17.000 


2 


1.937 


-1.076 


1.772 


3 


0.862 


190.070 


1.000 - 


4 


190.932 


145.121 


1.329 X 10 ^ 



Table 1 



FNR Equation ( 2 . 47 ) 


Iteration 


X 


Ax 


f(x) 


1 


3.000 


-0.667 


17.000 


2 


2.333 


-0.444 


4.160 


3 


1.889 


-0.296 


1.625 


4 


1.593 


-0.197 


1.124 


5 


1.395 


-0.132 


1.024 


6 


1.263 


-0.088 


1.005 


7 


1.175 


-0.058 


1.001 


8 


1.117 


-0.039 


1.000 


9 


1.078 


-0.026 


1.000 



Table 2 



Clearly the MNR equation causes x to diverge after an initial 
period of convergence while the FNR equation causes x to 
approach the minimum. 

6 . Computation Experience With the FNR and MNR Methods 
The performance of the MNR method in the preceding 
scalar example is typical of the performance observed by the 
author in large problems. However, the FNR equation is also 



\ 
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not uniformly effective when used exclusively in large 
problems. Fortunately, the areas of effectiveness of the 
two methods are complementary. 

In order to discuss the effectiveness of the two 
methods it is convenient to define two areas in the mlnlmiza 
tlon process. Initial behavior refers to the behavior of J 

a 

during the first two or three iterations in a sub-problem. 

Terminal behavior refers to the behavior of J after the 

a 

first two or three iterations within the same sub-problem. 
The following behavior has been observed. 

a. Initial behavior: The MNR equation outperforms 

the FNR equation in this area. The ability of the FNR 
equation to minimize J. is very sensitive to the starting 

a 

value of the unknowns (c). With the values of c far 
removed from the optimum, the FNR equation generally causes 
J. to Increase rapidly and diverge from the minimum. The 

a 

MNR equation on the other hand is relatively insensitive to 
the starting c and can usually be counted on to move J„ 
toward the minimum for at least one or two Iterations. 

b. Terminal behavior: As the minimum is approached 

the MNR equation produces the behavior shown in Figure 1. 

The FNR equation, however, generally becomes extremely 
effective in rapidly finding the minimum. 

7 . A Combination FNR-MNR Minimization Method 

The obvious approach suggested by the previous obser 
vatlons is to devise an algorithm which minimizes by the MNR 
equation initially in a given sub-problem and switches to 






\ 
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the FNR equation at some appropriate point In the Iteration 
process. Such an algorithm Is presented In Section IV. 
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III. AN INEQUALITY CONSTRAINT PENALTY FUNCTIONAL 



In this section a new penalty functional Is Introduced 
for state and control Inequality constraints of the form 

X 1 X (t) < X. , 1=1,2, ... ,I < n , t e [t ,T] , (3.1) 

Xl 1 Ijyj s - o 

u. 1 u.(t) < u , j=l,2,...,I < £ , t e [t^,T] . (3.2) 

Jl J Jjyi ^ ° 

All state and control Inequality constraints encountered In 
the problems solved herein are of this type. The difficul- 
ties encountered with existing penalty methods which led to 
the use of a new functional are related. The new penalty 
functional has performed well in computation and is used 
exclusively in the solution of the problems presented. 
Additionally, several desirable theoretical properties of 
the proposed penalty functional are presented. 

A. INTERIOR PENALTY METHODS 
1 . Past Research 

In Ref. 10 Jones and McCormick present a number of 
theoretical results concerning interior penalty functionals 
of the Flacco-McCormick type [Ref. 9] in conjunction with 
the epsilon method. If, for example, a state or control, 
denoted by y(t) for generality, is constrained by 

y(t) < Y , t e [t^,T] , (3.3) 
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a Fiacco-McCormick penalty functional [Ref. 9] of the form 



V '3.>o 

^O 

is added to the augmented performance measure. The behavior 
of the integrand of expression (3*^) for a fixed time 
teCt .T] as the positive weighting factor r approaches 0 
is shown in Figure 2. 




Figure 2 

Fiacco-McCormick penalty function vs. constrained variable 

for a fixed time 



If penalty functionals of this type are added to the perform- 
ance measure, it can be shown [Ref. 10] that as r approaches 
0 and e approaches 0, the epsilon method yields Pontryagin's 
minimum principle. The development parallels and augments 
Balakrlshnan* s work [Refs. 13 and l4] without inequality 
constraints. No computational results, however, are 
presented. 

2 . Computational Experienc e 

A simple problem Involving one state variable and 
one constrained control was attempted using the epsilon 
method and a penalty functional of the form given by 
equation The optimal control was on the constraint 

boundary. The algorithm was unable to solve the problem 
by either the FNR or MNR method from a variety of starting 
points. Once the control penetrated the constraint boundary 
for a finite time interval, the algorithm failed on the next 
Iteration. The value of r required to keep the control 
admissible for all t e [t^,T] throughout the iteration 
process was large, resulting in the augmented performance 
measure being dominated by the Placco-McCormick penalty 
term. As a result, the optimal solution [u*(t,e,r)] to the 
augmented problem could not be made to approach the optimal 
solution [u*(t)]. 

B. A NEW PENALTY FUNCTIONAL: COMPUTATIONAL PROPERTIES 

1 . The Form of the New Penalty Functional 

Consider a control or state y(t) which is subject to 
a constraint of the form 
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y(t) e [yj^,yj^] 



3 



(3.5) 



t e Ct^,T] 



A penalty functional of the form 

T 
o 

where is a positive integer is added to the augmented 
performance measure. The effect of Kp can be seen from 
Figure 3 which shows the integrand of equation (3.6) as a 
function of y(t) for a fixed time t e [t^,T]. 



Jp(y,r,Kp) = 



/ 



2y(t) - y^ -- yL ~|2K 



dt 



(3.6) 




Figure 3 

New penalty function vs. constrained variable for a fixed time 
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A functional of the form given by equation (3.6) is added 

to the augmented performance measure for each inequality 

constraint of the form given by equations (3.1) and (3'.2). 

For I control constraints and I state constraints of this 
^ s 

form the total augmented performance measure for the epsilon 
method written for time Invariant problems with t^ = 0 is 




f(x,u) 1 1^ dt 



(3.7) 




lM3,j2Kp| 



dt 



= J 




+ r J 

P 



(3.8) 



The "two-sided" feature of the penalty functional makes it 
especially suited to constraints of the form given by 
equations (3.1) and (3.2). In effect, two inequality 
constraints are Included in one penalty term. 

2 . Computational Strategy 

The power Kp is increased gradually in numerical 

computation in the same manner as e is decreased. Thus, 

increasingly refined boundaries to the admissible region 

are provided. Both e and Kp are held constant within a 

sub-problem and are altered between sub-problems. The 

weighting factor r, which is required to provide an overall 

weighting among J, J , and J , is held constant throughout 

s p 



the entire problem. 



Computational results with this penalty functional 
have been excellent. Up to eight inequality constraints have 
been treated successfully in one problem. 

C. THEORETICAL PROPERTIES OP THE NEW PENALTY FUNCTIONAL 
1. Introduction 

Three desirable properties of the new penalty 
functional are presented here. These properties and their 
importance are discussed followed by a proof of each 
property. 

a. Penalty functionals of the form of equation (3*6) 
are convex on R^ where c e R^ is defined by 

•s* O' 

= Ca]^,a2, . .. ,aj^^,y(t^),y(T)] (3.9) 

and 

y(T)-y(t ) M mir(t-t ) 

y(t) = y(y + — — <3. 10) 

o o 

It is desirable that the augmented performance measure be 
convex in the unknowns of the minimization to Insure that a 
global minimum is attained. If it can be shown that the 
Inequality constraint penalty functionals (3*6) are convex, 
then the addition of any number of these functionals does 
not destroy a convexity condition which exists without these 
terms, because the sum of convex functionals is also convex. 
Indeed, the addition of terms of the type given in equation 
(3.6) may create a convexity condition where one does not 
exist without the terms. 
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b. If for a fixed c given by equation (3*9)> the 

~y 

expansion (3.10) of a constrained state or control is inad- 
missible by a finite amount e (e >0) at t^eCt ,T), its 

p p o 

associated penalty functional (3.6) is unbounded as Kp 
This means that as Kp the contribution of a penalty 

term (3.6) to the augmented performance measure for an 
inadmissible state or control becomes very large. Therefore 
if J„ is being minimized under the condition of ever Increas 
Ing Kp, the constrained states or controls must at least 
approach admissibility. 

c. If for a fixed c given by equation (3.9)> the 

~y 

expansion (3.10) of a constrained state or control lies 
completely within the admissible region, its associated 
penalty functional (3.6) has limit zero as Kp The 

significance of this result is that penalty terms of the 
form given by equation (3.6) will add less and less to the 
augmented performance measure as Kp is increased for states 
and controls that are completely admissible . 

2. Convexity 

Property a discussed above is shown here. The 
theorem to be proved follows. 

Theorem 1 . If a constrained state or control y(t) 
is bounded for t e [t^,T] and is given by 

y(T)-y(t ) M mrr(t-t ) 

y(t) = y(t ) + — ^ (t-t ) + sin (3.10) 

^ o nFl o 



where c 

~y 



£ is defined by 
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c = [a^, aj^, y(t^), y(T)] 0 , (3-9) 



then the penalty functional 



/'[■ 



.Tr2y(t)-y -y n 
Jp(y.Kp) = r ' I « ^ 



^M-^L 



2K 



dt 



(3.6) 



,n 



is convex on R . The constants and define the 
admissible region for y(t) and r is a constant. 



Proof . Consider the case of K =1. Equation (3.6) 

P 

becomes 



J.^(y,K^) = — — J^l2yit)-iy^+y^)T dt. (3.11) 



^ ^ ( V “V ) ^ t 

vyM yL-» 



Let 



A ^L 

d 2 



(3.12) 



and 



r 

o 



Hr 






(3.13) 



Substituting equations (3.12) and (3.13) into equation (3.11) > 
we obtain 



37 



dt 



(3.14) 






T 

Cy(t) - d] 



2 



= r 



o 




- 2y(t)d + d^] dt . 



(3.15) 



Substituting the expansion (3.10) into equation (3.15) > we 
obtain 






y(to) + 



y(T)-y(t^) 

f=t 

o 



M 

“-‘o’ ‘ k ^ 



sin 



imr(t-t^) 

o 



( 3 . 16 ) 



r y(T)-y(t ) 

L o 



M 

L 



sin 



iniT(t-t ) ■ 



d + d' 






dt. 



Equation (3.16) may be written as a quadratic functional 
of the form 



J 

P 



o 



T 

f [i-<£y’Sl^^^£y> + (3^17) 



O 



where a( t ) e is 




= -2d £ 



sin 



Tr(t-t^) 
~ T-t 

o 



sin2Tr(t-t^) 



sinMTr(t-t ) t-t t-t n 

o n 9.) 2. 

T-t^ T-tJ »T-t^ J * 

(3.18) 






38 



Qj^(t) is the outer product given by 



Q (t) = aMi^ 



( 3 . 19 ) 



eind 



3 = . 



( 3 . 20 ) 



At an arbitrary fixed time t*e[t^,T], the integrand of 
equation (3.1^) is 



Cy(t^f)-^d]^= I" < £y ^ 



( 3 . 21 ) 



The first term on the right side of equation (3.21) is 
1 , , r y(T)-y(t ) M mTT(t^-t ) -|2 

5<2y.«l(t.)Cy> = [y(t^) . ^ ^ • 

( 3 . 22 ) 

Since the terms in the finite expansion of y(tjj) given in 
equation (3.10) are linearly independent and analytic, y(t^) 
is different from zero almost everywhere and 



y‘^(t«) > 0 , t« e [t^,T], c,^ 7 ^ 0 



~y 



( 3 . 23 ) 



Thus , 






( 3 . 2 ^^) 
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and is, therefore, positive semi-definite, at least, 

and positive definite almost everywhere for any c 0. 

~y 

Applying Theorem 4.5 of Reference 17 (p. 27)^, the function 
given in equation (3.21) is convex on at t = t^. 

Next, consider the case where is any positive 

integer. In Appendix E the following theorem is proved: 
if f(x) is convex on R^ where x e R^ and f(x) ^ 0, then 
f^(x) is convex on R*^ where K is any positive integer. 

Since 



[yCt,^) - d]^ > 0, (3.25) 

is convex, it follows immediately that 

2K 

Cy(t^) - d] P (3.26) 

is convex on R^ at a fixed time t« e [t ,T]. Since y(t) is 

s O 

bounded by assumption for t e [t^,T], it follows that for 



Let f be a twice continuously differentiable real- 
valued function on an open convex set c in R . Then f is 
convex on c if and only if its Hessian matrix 



Qx = (QijCx)), 



/ 'V _ 3^f 



(C 



1’ 



* ^n^ 



is positive semi-definite for every x e c. A quadratic function 



f(x) = % <x,Qx> + <x,a> + a 

where Q is a symmetric n x n matrix, is convex on R^ if and 
only if Q is positive semi-definite. 



any finite positive integer K , the expression (3.26) is 

P 

bounded for tj^ e [t^,T]. By Theorem i|, (p. 536) of 
Reference 18 ^ 




(3.27) 



is convex on R^. Hence equation (3.6) is convex on R^ for 

all finite positive integer values of . 

3 . Behavior of the New Penalty Functional for an 
Inadmissible Constrained State or Control 

Property b is shown below. 



Theorem 2 . Assume y(t) is bounded and given by 

equation (3.10) for t e [t ,T] where c 7 ^ 0, as defined by 

o -y -- 

equation (3.9). If lor a given y(t) 



> yj/i (3.28) 

or 

y(tjf) 1 yL - £p (3.29) 



^Let 

Ip(x) = 

Let T be of finite measure. Let f(t,x) be a finite convex 
function of x for each t and a bounded measurable function 
of t for each^x. Then is a well-defined finite convex 
function on l'°(T) which Is everyv;here continuous with 
respect to th§ uniform norm. 



,x(t)) dt. 



at some time t^^ e (t^,T), 



where e > 0 
P 



then 



Limit J (y,K ) = « ( 3 . 30 ) 

K ^ P P 
P 



where 



J (y,K ) 

p * p 



r 




2y(t)-y^-yL-| 

yM-^L 




dt . 



(3.6) 



In order to prove this theorem the following lemma is 
required. 



Lemma 1. Assume y(t) is bounded by 



-00 < < y(t) < M2 < ~ • (3. 31) 

and is given by equation (3.10) for t e [t^,T] where 
c 7^ 0 as defined by equation (3.9). Then, if 

y(t*) 1 yM ^p (3.28) 

at some time t^ e (t^,T) for any > 0, there exists a 
6 > 0 such that 

y(t) > y^^ + (3. 32) 

for the finite time interval 

'S£t£tjf+ 6. 

l\2 



(3.33) 



Similarly, if 



y(t^) < yL - Ep (3.29) 

at some time t- e (t^,T) for any e >0, there exists a 

o p 

6 > 0 such that 



y(t) < y^ - ^ (3.3^) 

for the finite time interval 

t^f - 6 < t < tjj + 6. (3.33) 

Proof of the lemma . First, it is necessary to show 
that the coefficients in equation (3.10) are bounded; that 
is 



|a^| < m = 1,2, . . . ,M 



(3.35) 



where > 0. To this end consider 



/ „ iUp y(T)-y(t ) M mn(t-t )-i2 

[y(t)] dt = jy(t^)-i — 5:5- — — J 



dt 



(3.36) 



= < c ,Q^c > 



(3.37) 



^<3 



I \ 



/r! 



where c„ is given by the definition (3.9). Qo is a 

2(T-t ) 
o 

ir 

_ 2(T-t^) 

2i 

« 

« 

” Mir 
3 
■■'3 

( 3 . 38 ) 

where the terms with i are positive if M is odd and negative 
if M is even. Since y(t) is bounded by assumption and the 
expansion (3.10) is the sum of M + 2 linearly Independent 
terms, we have 

T 

0 < f [y(t)]^dt £ ( 3 . 39 ) 

t*' 

o 

where > 0. Using equation (3*37) and inequality (3*39), 
we obtain 



~y 






symmetric matrix given by 









?2 



T-t. 

~T 

0 



b 



2(T-t ) 
o 



TT 



2(T-t^) 

TT 



0 ... 0 

T-t n 

o . .0 



2(T-t^) 

2tt 



2ti 



T-t 



o 



2(T-t^) 

Mi 



_2(T-t^) .+2(T-t^) 



Mti 



2(T-t ) 
o 

TT 

2(T-t^) 

2i 



2(T-t^) 

Mi 

T-"o 

3 

T-t 



0 < <3y>S2Sy> ^ 






\ 



( 3 .^ 0 ) 



^2 is, therefore, positive definite. Using Theorem 2.5 
of Reference 15 (p. 52)^, we have 

<2y’S2Sy> - ^ Il2yl 1^ (3.^1) 

where ^ > 0 is the smallest eigenvalue of and 
ll£y|| = yj<^ 3y »Sy> • Therefore, using the inequality 
(3-^0), we obtain 






(3.^2) 



Since a_, m = 1,2 
m ’ 



,M is a subset of 



c 



~y 



> 



it follows that 



l^ml - ^3’ m=l,2, .. . ,M 



( 3 . 35 ) 



where > 0. 

Now consider the derivative of equation (3.10) 



which is 



y(t) = 



y(T)-y(t^) 



O 



ITITT 



T-t 



O 



M 

„?T^m T-t 
m=l o 



cos 



mir(t-t^) 



T-t 



o 



( 3 .^ 3 ) 



^Let Q = (^ii) ^ symmetric n x n matrix. Then Q is 

positive~definIte if and only if there is a k > 0 such that 

< v,Qv> > k| |v| 

for all V in R^, where | |v| | = •»/< v, v > 

norm of v. 



is the Euclidean 



Taking the absolute value of both sides of equation (3.43) 
and applying inequality laws , we obtain 



. y(T)-y(t ) 

|y(t)| < I ^ 



, , I. mn M(t-t ) 

L?-, ■” — 

m=l o o 



(3.44) 



which further simplifies to 



y(T)-y(t ) M 

|y(t)l 1 I T-t ° I ^ 2 1^1 

o m=l o 



(3.45) 



Applying the Inequality (3.35) > we obtain 



.. y(T)-y(t^) M,tt M 

IHt)i i I — T-t I ^ T^r- 



'o m=l 



< 

■" 5 



(3.46) 

(3.47) 



for all t e (t ,T) where Mr- > 0. The first part of the 
o 5 

lemma as expressed by the inequality (3.32) will now be 
shown. Let 



6 ^ -£ — 



(3.48) 



as shown in Figure 4. Consider 






• « i • 




M 



1 





i 




Applying inequality (3*^7), we have 



y(t) ^ y(t*) - Mf. Itjf - 1 1 . 



(3.50) 



Consider t in the interval t^^ - 6 £ t £ tj^ + 6. In this 
interval 6 satisfies 



5 i |t* - t 



(3.51) 



Applying inequality (3.50), we have 



y(t) ^ y(tjf) - 6 



(3.52) 






Using the definition (3.^8), we obtain 



y(t) ^ y(t*) 

Applying inequality (3.28), we 

1 yjvj + 



e 

- . 

2 * 

obtain 



. e 




or 







( 3 . 53 ) 



( 3 . 5^0 



(3.32) 



thus proving the first portion of Lemma 1. The second portion 
of the lemma as given by inequality (3.3^) can be proved in 
a similar manner. It is possible at this point to return 
to the proof of Theorem 2. 



Proof of Theorem 2 . If inequality ( 3 . 28 ) applies, 
then by Lemma 1 inequality (3.32) is true. Considering the 
integrand of equation (3.6), since y^^j > y^^, it follows that 



• 2y(t)-y^-yj.- 
^M"^L - 



2K 



> 



y^yj yj;^ ' 

I ^m”^l 




( 3 . 55 ) 



for tjj-6_<t£tjf+6. Further simplification yields 



[■2y(t:)-yf,-yL 


r'^p > 1 


rs^^M-^L 1 


2Kp 




- 1 


L yr/i-^L J 


y 




> r 


1 + ^P-- 






-1 


^M-^L 


J 


Since the integrand 


of equation (3.6) 


is nonnegative 


t e (t^,T) and r > 0 


, it follows that 




0 


j"2y (t )- 




dt 


1 ^M- 


-^L J 



(3.56 



(3.57) 



(3.6) 



> r 



T ,, , 

tj-s L j 



'M -^L 



(3.58) 



Applying inequality (3.57) > we have 



t^+6 



Jp(y,Kp) > r 



./. [ 



1 + 



yM-^LJ 



2K 



dt . 



(3.59) 



The Integration of inequality (3.59) yields 



Jp(y 



.Kp) 1 r [ 



e„ h2K 



1 + 






26 . 



( 3 . 60 ) 



Since yj^ > y^^ and e > 0, it follows that 



Limit 
P 



K_ •.» 1^ %% 



2K 



P = « 



( 3 . 61 ) 



^9 



and in view of inequality (3.60) 



Jp(y.Kp) = ” • 



(3.62) 



If inequality (3-29) applies, then by Lemma 1, 
inequality (3-3^) is true. Since yj^^ > y^^, it follows that 



[ 



2y (t)-y^- y^ 



14 



2(yL- -#)-yM-yLi 



^M-^L 



(3.63) 



for tj(-6j<t£tjf + 6. Simplifying, we obtain 



1 

ct 

f 

1 


1 - r 1 i 1 


L yM-^L 


J - L yM-yL J 



(3.64) 



Squaring both sides of inequality (3.64), we have 



r ^y(t)-yM-yL f , r, , 

I- ^iC^L J - L J 



(3.65) 



Raising inequality (3.65) to the th power, we obtain 



2y(t)-y^-y^- 

^M“^L 



2K 



r e l2K 

> I 1 + — f 



(3.57) 



which is identical to Inequality (3-57). The remainder of 
the proof follows inequality (3.57) to equation (3.62) 
exactly. The theorem is proved for the open interval t e (t^,T) 
The extension of the theorem to cover the closed Interval 
t e [t^,T] is not difficult. 
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^ • Behavior of the New Penalty Functional for an 
Admissible Constrained State or Control 

Property c above Is shown below. 

Theorem 3 . Assxame y(t) Is given by equation (3.10) 
for t e [t^,T]. If for a given y(t) 



1 1 - =q 



( 3 . 66 ) 



for all t e [t^,T] where > 0, then 



Limit / TT \ r\ 
^ J(y,K)=0 

Kp 00 p ' » p " 



(3.67) 



Proof. Prom Inequality (3.66) It can be seen that 






( 3 . 68 ) 



since yj^ > y^^. Inequality (3.68) may be rarranged to the form 



1 



2e, 



- 



(3.69) 



The Inequality (3.69) will become useful shortly. 

Starting with Inequality (3.66), multiplying through 
by 2, and subtracting yjyj and yj^, we obtain 



- 2y(t)-yj^-yL 1 2(yi^-eq)-yiy[-yL • (3.70) 
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Dividing through by the positive quantity y^^ - yj^, we obtain 



^<yL^e^q)-yM-yL ^ ^ 2(yM-eq)-yn-yL 



(3.71) 



This inequality can be reduced to 
2e 



- 1 



— ^ < r 

yw-^L - L 



2y(t)-yj^-yLl 

yM-^L 



2s 

< 1 3L 



yM-yL 



(3.72) 



In view of inequality (3.69), inequality (3.72) may be 
rewritten as 



i[ 



2 y(t)-y„-yj^ 

^M-^L 



]U 



1 ^ 



yjvi 



(3.73) 



Raising inequality (3.73) to the 2Kp power, we have 



|C 2y(t)-yM-yL ' 

I L yM-yL 



2K 






■|2K 
^1 P 



yjyi yL-l 



(3.7^) 



Observing that 

■2y(t)-yj^-yLn 2K 



ir 



yM-yL 



p _ |- 2y(t)-yj^-yL -i2Kp 

L yM-yL J 



(3.75) 



we have from inequality (3.7^) 



r2y(t)-yj^-yLi 


'"p < 


( 

a 

CM 

1 

r~ 


L yM-^L J 


v» 


yM-yL -1 



2K 



(3.76) 



for all t G [t^,T]. In view of inequality (3.76), it is 
easily seen that 
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2y(t)-yj^-yj^-|2K 



dt 






Performing the integration, we have 



■2y(t)-yj^-yj. 



r"r2y(t)-y 

J L 



M 



2K 



dt < 



2e 



1 - 



q 






2K 



(T-t^) .( 3 . 78 ) 



In view of Inequality (3* 69 ), and the fact that yj^ > yj^, 
it follows that 



Limit r, 

r 



2e. 







0 . 



Observing inequality (3* 78), we obtain 



(3.79) 



Limit 

Kp-«> 



,/ 



V 2y(t)-y^-yj.- 
l- yM”^L 



2K 



dt 



0 



( 3 . 80 ) 



The theorem is proved. 
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IV. THE ALGORITHM 



This section describes the algorithm used for minimizing 
the augmented performance measure . 

A. GENERAL MINIMIZATION STRATEGY 

1 . Sequence of Unconstrained Sub-problems 

A sequence of unconstrained sub-problems is solved 
by the algorithm. In each sub-problem the algorithm mini- 
mizes the augmented performance measure for given values 
of the weighting factors (e and r) and the inequality con- 
straint penalty term power (K ) . After an appropriate 
stopping criterion is satisfied, e is reduced, is 
increased, and a new sub-problem is commenced using the 
optimal solution to the last sub-problem as a first guess. 
This procedure is repeated until enough sub-problems are 
completed to meet a second stopping criterion. 

The algorithm is programmed to do one sub-problem 
on each computer run. The results are stored on an external 
storage device between computer runs and are retrieved at 
the commencement of the next run (new sub-problem) . 

2. Minimization Strategy 

The algorithm minimizes by either the FNR or MNR 
method. The user must decide which method to use on each 
iteration. This is a matter of experimentation, especially 
for the first two or three sub-problems. An effective pro- 
cedure is to run the sub-problem once using the MNR equation 
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throughout and once using the FNR equation throughout. From 
these results an effective minimization strategy can generally 
be deduced for the sub-problem. Occasionally further exper- 
imentation is required. This experimentation points out 
the advantages of using separate computer runs for each 
sub-problem. Once a sub-problem is completed and the results 
stored, the computation does not have to be redone each 
time an experimental run is made in the next sub-problem. 

B. COMJ/IENCING THE PROBLEM 
1 . Initial Decisions 

Three interrelated decisions must be made to begin 
a problem. First, the number of time points K must be chosen. 
Second, the number of coefficients M for each state and 
control expansion must be chosen. The same number of coeffi- 
cients is used for all expansions in a given problem in 
this dissertation, but this is not a requirement. From a 
theoretical standpoint it is desirable to use a large number 
of coefficients and time points to insure that an adequate 
approximation of the optimal control and state trajectory 
is obtained, but practically, computer time and storage 
requirements limit the number of each. The computational 
penalty for using a large number of coefficients is the 
more severe of the two as the number of equations in (2.36) 
and ( 2 . 37 ) which must be solved is equal to the total number 
of coefficients plus the number of free end conditions. The 
solution of equation (- 2 . 36 ) or (2.37) represents a considerable 
portion of the overall computer time. 
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The third decision involves the initial values of 

e, r, and K^. The weighting factor r for the inequality 

constraint penalty terms is held constant throughout the 

entire problem. A satisfactory value used in all problems 

in this dissertation for all inequality constraint penalty 

terms is r = 100. An initial value of K which has worked 

P 

well in all problems is = 4. Larger values of gener- 
ally cause computer overflow in the first sub-problem. 

With these values chosen there exists a region of e’s for 
which the first sub-problem will respond to an appropriate 
minimization strategy. This acceptable range of starting 
e's is different for each problem but is in the range 

10“^ < e < 10“^ 



for all problems solved herein. Numerical experimentation 
is the only method available to determine an acceptable 
starting e. There is no theoretical requirement to use 
the same value of e for each state equation equality con- 
straint term in the augmented performance measure or the 
same in each inequality constraint term, but the use of 
different e's and K^'s has never been required. 

2. Initial Guess for the Unknowns 

Once the above three Initial decisions are made, an 
initial guess for the vector of unknowns c is required. The 
vector c includes all coefficients and free end conditions. 



X ^ 
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All coefficients are set equal to zero initially unless 
there is good reason to make a different choice. 

C. ITERATION 

1 . Required Vectors and Matrices 

The states and controls are calculated at each time 
Increment by evaluating the functional expansions. The w 
vector defined in (2.15) is calculated using these states 
and controls. Next, the gradient matrix (2.3^), the aug- 
mented performance measure (2.l6), the symmetric matrix 

T T 

(V^w)^, and the vector (v w)„ w are calculated. 

At this point the algorithm begins the iteration 

process with either the MNR or the FNR method depending on 

the value of a flag set by the user (the method selected 

is based on the iteration number being performed). If the 

MNR method is to be used, equation (2.37) is formed. If 

the FNR method is called for, the three-dimensional array 

2 T 

(2.35) is calculated and the symmetric matrix (V w) w 

c ^ c ^ c 

is formed. It is prohibitive to store the entire three- 
dimensional array, but a feasible alternative is to multiply 
each matrix in this array by w as the matrix is calculated 
and store the resulting coliiinn vector. Once a matrix in 
the three-dimensional array is multiplied by w, it is no 
longer required by the algorithm. The next matrix in the 
array is calculated and stored in the same storage locations 
used by the first matrix. Only the symmetric matrix 

P m 

(V w) w need be stored. The total increase in storage 
c ~ c ~c 
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requirements of the FNR method over the MNR method using 
this computation technique is less than 10 percent in the 
problems solved herein. It is also imperative in terms of 
computation time to take full advantage of the symmetry of 
the matrix (7 w) w . Due to this symmetry it is necessary 

C ^ o 

to calculate only one column of the first matrix in the 
array, two columns of the second matrix, and n columns of 
the n^ matrix. By taking advantage of the symmetry the 
average time for each FNR iteration is approximately twice 
the time for each MNR iteration. 

2 . Solving the Linear System 

At this point equation (2.36) or (2.37) is formed 
and must be solved for Ac. This is a linear system of the 
form 



A X = b (i|.l) 

and is solved in the algorithm by one of three methods 
available to the user. They are: 

a. Gauss elimination with Improvement by residuals 
using total pivoting, 

b. Gauss elimination with improvment by residuals 
using main diagonal pivoting and a computation technique 
which capitalizes on the symmetry of A, and 

c. Gauss-Seidel iteration. 

In the problems solved herein the number of unknowns 
varied from 37 to 7^. In spite of the large number of 
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unknowns Involved, the elimination methods required less 
computation time to solve the linear system than the Gauss- 
Seidel iteration method. It was observed for the problems 
solved that total pivoting was not required in the elimina- 
tion method. Method b, therefore, was the most economical 
and effective method for solving the linear system and was 
used in all problems. Method c is retained in the event 
that the algorithm is used to solve problems with a larger 
number of unknowns . 

In each solution of equation (M.l) one Improvement 
is made using residuals. That is, after equation (4.1) is 
solved. 



A X - b = r (4.2) 

is formed. The system 

A y = r (4.3) 

is solved and the resulting y is subtracted from x to form 
the final solution to equation (4.1). 

3. Interpolation 

Tabular functions of two independent variables are 
used extensively in the problems to represent aircraft and 
missile parameters accurately. Parabolic interpolation is 
used to obtain the functional values in these tables and 
the required first and second partial derivatives. The 
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derivation of the necessary difference equations for para- 
bolic interpolation in two Independent variables is presented 
in Appendix C, 

Excerpts from the tabular data used in the problems 
is presented in Appendix B along with graphical representa- 
tions of the data. The data represents typical supersonic 
aircraft and missile performance parameters and has been 
obtained from several sources. Considerable effort was 
expended to smooth the data before the tables were constructed 
since finite difference methods were used not only for 
functional values but also for first and second partial 
derivatives . 

^ . Stopping Criteria 

Once the linear system is solved, a new c vector is 
calculated from 



c^'^^ = + Ac^ . (H.H) 

At this point a stopping criterion is tested. If 

|J ^ < STOPl , (4.5) 

Si Si ^ 

the sub-problem is finished. Otherwise the iteration process 
is continued. At the completion of the sub-problem the 
results are stored off line. The computer run is complete. 

To begin a new sub-problem a new computer run is 
initiated, recalling the results stored from the last 
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sub-problem. Epsilon is decreased, K is increased, and 
the minimization strategy is altered by the user as required. 
Typically, e is divided by a factor of between two and ten 
and K is increased by two or four, That is, 



K = K i + 



2 

or 

i\ 



(^. 6 ) 



More ambitious policies usually result in failure of the 
algorithm. 

Sub-problems are solved until a second stopping 
criterion is satisfied, Several criteria are possible to 
end the problem, A method used successfully Involves 
observing 




and stopping when this sum, which represents the' penalty 
terms due to the equality and inequality constraints without 
weighting factors, ceases to decrease significantly between 
sub-problems . 

5 , Flow Chart 

A flow chart of the algorithm is given in Figure 5. 

6 , Integration 

At the completion of the last sub-problem a check 
on the degree. of satisfaction of the state equations is 
obtained by comparing the state expansions with the state 
trajectory obtained by integrating the state equations with 
the control expansions. 
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Figure 5 

Algorithm Flov; Chart 



\ 



62 




V. A MISSILE INTERCEPT PROBLEM 



In this section a short-range missile intercept problem 
is solved. An air-to-air missile launched from an attacking 
airplane is to intercept a constant-velocity target in 
minimum time. The missile is restricted to move in a plane. 
The orientation of this plane is defined in three- 
dimensional space as the plane containing the position of 
the missile at launch, the position of the target at launch, 
and the velocity vector of the target. The assumptions 
applied to the problem, the coorinate systems used, the 
nomenclature, and the derivation of the equations of motion 
are presented in Appendix A. 



A. PROBLEM FORMULATION 
1. State Equations 

The state equations derived in Appendix A are 



• gTh gpSa 2 gpSa 2 SW 

M = Th cosa - 2W~~ ^ “ W~~ ^ ~ ivT 

e e e e 

(5.1) 



0 = 



gThjj^ Th sina 



aW, 



M 



gpSa gpSa 

2W f^C^sina + 



MC„cosa 

N 



gW cos 9 
^ c 



aW, 



M 

(5.2) 



aMcos0 


- aM^cosy 




(5.3) 


aMsin0 


- aM^siny 


• 


(5.^^) 
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The states are Mach number M, flight path angle 6 , relative 
range X, and relative cross-range Y as defined in Figure ^ 0 . 
The control is angle of attack «. 

The normalized missile thrust Th is considered 
constant until missile engine burnout tg and zero thereafter. 
That is. 



Th = 1 , t e [ 0 ,tg] 

Th = 0 , t e (tg,T] . 



Other parameters considered constant are 



g 

Th 

W 

€ 

s 

t 

a 

P 



M 



B 



32.1725 


ft . / sec , ^ 


3500 


lbs. 


200 


lbs . 


0.35 


ft.^ 


7 


sec . 


1077.8 


ft . /sec. 


0.001756 


slugs/ft . 



( 5 . 5 ) 

( 5 . 6 ) 



( 5 . 7 ) 



The values of the constants in equations ( 5 . 7 ) are based on 
an engagement at 10,000 feet altitude. It is convenient to 
group these constants as 



C 



1 



STh„ 




( 5 . 8 ) 



gpSa 



(5.9) 



C 



2 




C 



3 




( 5 . 10 ) 



In addition, the computer solution is aided by normalizing 
the states X and Y by defining 



X 




(5.11) 



y = 



A Y 



( 5 . 12 ) 



where X and V are assigned nominal values of 10,000 feet. 
With these adjustments the state equations are 



M = C^Thcosa - C 2 M C^cosa - Cj^sina - C^W^sinS (5.13) 

» rp}, ci-fnrY W COS0 

® g - ” - C2MC^slna + C2MCj^cosa - (5.l4) 

X = ^ M cose - Y I^ipCosY ( 5 . 15 ) 

y = ^ M sine - y M^siny . (5.l6) 

2. Tabular Functions 



The axial and normal force coefficients and Cj^ 
are given in Appendix B as tabular functions of Mach number 
and angle of attack. This data is based on a typical air-to- 
air missile. 
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3 . Performance Measure 

The performance measure for this problem is 

J = / dt . (5.17) 

0 *' 

^ . Inequality Constraints 

Pour inequality constraints are required. The angle 
of attack must satisfy 



Prom structural considerations the load factor must satisfy 



5 . End Conditions 

In order to describe the initial conditions for the 
problem it is necessary first to pose the problem in the 
(X,Y,^) coordinate system shown in Plgures 4l and 42 of 
Appendix A. The problem chosen for presentation in this 
section involves a target in a shallow climb at short range 
with a slight altitude disadvantage crossing the attacker's 
flight path extension at 90°; l.e. 




( 5 . 18 ) 



-50 < I (0M) < 50 . 

o 



(5.19); 



= 15000 ft. 
hy = -3000 ft. 

3^ = 90° 



( 5 . 20 ) 

( 5 . 21 ) 

( 5 . 22 ) 

(5.23) 
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Following the procedure outlined in Appendix A, the 
remaining unknown parameters are 



= 0.8 (5.24) 

Y = 42.35° (5.25) 

W = 51.526 lbs. (5.26) 



The initial conditions as computed by this procedure are 



M(0) = 0.8 (5.27) 

0(0) = -49.573° ( 5 . 28 ) 

x(0) = 0 ( 5 . 29 ) 

y(0) = 0 . (5.30) 

The final conditions as computed by this procedure are 

x(T) = 0.9920 (5.31) 

y(T) = -1.1645 . (5.32) 



Note that x and y are the normalized relative range and 
relative cross-range, respectively, as defined by equations 
(5.11) and (5.12). The states X and Y in equations (5.11) 
and (5.12) are defined as the position of the missile in 
the (X,y) coordinate system shown in Figure 40. At t = 0 
the missile is at the origin of the (X,Y) system, hence 
equations (5.29) and (5.30) apply. At t = T the missile 
must intercept the target. Since the (X,Y) system becomes 
fixed with respect to the target at launch, x(T) and y(T) 
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are equal to the normalized coordinates of the target In the 
(X,Y) system at launch and are given by equations (A. 93) and 
(A. 95). 



B. THE EPSILON METHOD FORMULATION 

1 . The Augmented Performance Measure 

Using the penalty functionals described In Section 
III for Inequality constraints, the augmented performance 
measure Is 




where e and r are weighting factors and 6 Is a constant used 
to make minor adjustments In the admissible regions. In all 
problems 6 Is given a value of 1.03. This adjustment Is 
applied to all admissible regions of constrained states and 
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controls. The power is limited in computation to approxi- 
mately 30 before computer exponent overflow problems develop. 

With this value of K it is desirable to adjust all admls- 

P 

slble regions slightly as can be seen from Figure 6. 




Figure 6 

Adjusted admissible regions 
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The required elements of w are 



”k = [ 
'^K+k 

t 

^^2K+k 

'^SK+k 
"ilK+k ^ 

'^5K+k 

^6k+1 

where 



’k - Vhoosc.^ + 



I- CjW^slneJ , k = l,2,...,K 



(5. S'!) 



f«k - =1 



Thsina 



k 



^ =2’’k^A^=^"“k - °2”k°N, °°=“k 
k k k 



(5.35) 



w cose 



h c _£_ki r ^1^ k = 1 2 

3 M, J I G J * • • • * 



K 



= - I M^cose^ + I M^cosy] , k=l,2,...,K (5.36) 

' [^k “ 7 V^"^®k 7 [¥T^ > k=l,2,;. ..,K (5.37) 



r2a. "iK , 

' ^ ' 1,2, ...,K 



rae, M. “iK . 

^ [solrj - k=l,2,.... 



K 



(5.38) 



(5.39) 



= [(K-l)At] 






(5.!)0) 



\ = M(t) 



t=(k-l)At 



etc . 
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2 , Functional Expansions 



The state and control expansions written In discrete 



form are 



= Ml + (k-1) + £ a^^sln - - , k=l,2,...,K (5.^1) 

m=l 



0T,-0, M /-v 1 ^ 

®k " ®1 ~ Y-l ' ' > k=l,2,...,K. (5.^2) 

m=l 



^k = 



Xr-Xi 



M 



^1 K-l'" ' ^ . k=l,2,... 

m=l 



,K (5.^3 



^ i:gi (k-1) ^ f; a^sin • X=1.2.---.K (5.M) 

m=l 



M 



mu (k-1) 



“k ' “i ^ -fcr 

m=l 



=1,2, ...,K. (5.^15) 



3. Vector of Unknowns 



The elements of the vector c are 



( a^j^ , 9-2 j • • • j ®’]yi > ^2 * * * * *^ 



jyjjC^ J <^2 J • • • 



> ^2 > • • • > * * * ’ ^M’ ’ *^1 * *^K ’ ^ 



4. Partial Derivatives of the Tabular Function 



(5.i»6) 



The values of C. and C., are obtained from the tables 

90. 

by parabolic Interpolation along with the values of -rTyr- , il , 

o rl 9 ot ^ 



71 









90^ 9C^ 


— T~ * 


7~T~ 


9 9 


9 


'dVr 


da 


9M9a 


9M 9a 


2 








* 


The 


procedure 


Is outlined 


5. 


First 


Partial 


Derivatives 



2 

d‘^C 



3M 



N 

2“ » 



3^C 



N 



and 



9a‘ 



The first partial derivatives Indicated In equation 
(2.3^) are required. These partial derivatives are easily 
obtained from equations (5.3^) thru (5.^0) by taking the 
partlals of these expressions with respect to the vector c. 
The expressions are too numerous to Include. A typical term 
Is 



9w, 3W, 3M. 9 w 



(5.47) 



For 1 < k < K, Is given by equation (5.34) and the partial 
derivatives Indicated above are 



9mT 



' [t]”* 



(5.48) 






3C„ 

+ =2\ mr 



(5.49) 



“k ] m 



3M 

Ta 



k mtr • mTr(k-l) 

(K-l)At K-1 



(5.50) - 



m 
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(5.51) 



3M, 

t 

3a 



= sin 



m 



mu (k-1) 
K-1 



Notice that C. and C., are functions of M. as well as a, . 

6 . Second Partial Derivatives 

The second partial derivatives indicated in the three 
dimensional array (2.35) are also required. Again the expres 
sions are too numerous to list. A typical term for 1 < k < K 
is 



a "k . 9_ 

Lsa„J ' 



9w i 

Letting — — = R , we have 
3 



25- 9R. !!Ik + as_ !!lk 



The partial derivatives indicated above are 



(5.52) 



(5.53) 



3R _ 



3l^k 



3^ 

3M,. 



= 0 



(5.5^) 



■t 



3C, 



3^C, 






- costtj^ + 



3M. 



— cosa^^ 



k 



+ 2CoC.. sina,, + 
2 N|^ k 



9C„ 9^C„ 

2 "k 

2-k 9M7^ =i"“k " ^2«k TTT^ 



sina, 



3M, 






'Ssin !Sl(j£:dI 



K-1 



(5.55) 
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_ l-n ZttOc-1) 

(K-l)At K-1 



(5.56) 




(5.57) 



C. RESULTS 

The problem was solved twice: once using 8 coefficients 

for each expansion (problem A) and once using 12 coefficients 
for each expansion (problem B). In both cases K = 21 time 
points (20 time intervals) were used. The initial guess 
for the c vector was: ' 



1 . Problem A - 8 Coefficients for each Expansion 

Pour sub-problems were required to solve the problem. 
Table 3 gives the weighting factor values, optimization 
strategy, and computer time for each sub-problem. Table ^ 
gives the components of the minimum augmented performance 
measure for each sub-problem where 



all expansion coefficients = 0 



M(T) = 1.4 
e(T) = 0® 

a(0) = 20® 
a(T) = 0® 



(5.58) 



At = 7/20 sec. (T = 7 sec.) 




(5.59) 
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sub- 

problem 


c 


R 


K 

P 


I* 


optimization 

s,tJ?ategy** 


C.P.U. 

time 


1 


1.0 X 10"^ 


100 


4 


8 


MMMMMMMM 


2 ’28" 


2 


0.67 X 10"^ 


100 


6 


6 


FFFFFF 


3' 38 " 


3 


0.5 X 10"^ 


100 


8 


2 


FF 


115411 


1} 


0.5 X 10"^ 


100 


32 


1 


F 


1' 32" 



* number of iterations required 
** M - MNR method 
F - FNR method 



Table 3 

V/eighting factors, optimization strategy, and C.P.U, time for 
missile Intercept problem A. 



sub- 

problem 




J* 


* 


* 


1 


16.95 


7.538 


0.6853 X 10 "^ 


0.2554 X 10 "^ 


2 


14.52 


7.585 


0.4626 X 10 "^ 


0.1287 X 10 "^ 


3 


16.83 


7.585 


0.4624 X 10 "^ 


0.4619 X 10 "^^ 


4 


16.83 


7.585 


0.4624 X 10 "^ 


0.3691 X 10 "^^ 



Table 4 

Components of the minimum augmented performance measure for 
missile Intercept problem A. 
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Augmented performance measure 



Figure 7 is a plot of the augmented performance 



measure vs. iteration number. 




Figure 7 

Augmented performance measure vs, iteration 
number for missile intercept problem A 
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Iterations performed by the FNR method are Indicated by a 
solid line. Iterations performed by the MNR equations are 
indicated by a broken line. The figure shows several signi- 
ficant characteristics: 

a. the failure of the FNR method on the first itera- 
tion of the problem (the initial guess is too far from 
optimum to allow the FNR method to converge); 

b. the superior terminal convergence produced by 
the FNR method close to the minimum; 

c. the oscillatory results produced by the MNR 
method as the minimum is approached. After the commencement 
of sub-problem 2 as shown in Figure 7, the FNR method is 
used exclusively to the end of the problem. For comparison, 
sub-problem 2 is commenced with the MNR method and allowed 
to run to 23 iterations. At the beginning of sub-problem 3> 
J. increases slightly and never returns to the minimum value 

ci 

obtained in sub-problem 2. This is an indication that e has 

reached a point where smaller values have little Influence 

on the reduction of J . 

s 

Figure 8 is a plot of the performance measure vs. 
iteration number. The performance measure (final time) 
Increases with Iteration number. As the algorithm Iterates, 
At is being Increased to reduce (the term in the augmented 
performance measure reflecting the degree of non-satisfaction 
of the state equations) while insuring that constrained 
states and controls are kept within admissible bounds by 
reducing or holding down J . With small values of e and 
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Sub-problem number 




Figure 8 

Performance measure vs. iteration number 
for missile intercept problem A. 

large values of the end result of minimizing the augmented 
performance measure is a control history and state trajectory 
which minimizes the time to intercept while satisfying the 
state equations and inequality constraints; all to a 
reasonable degree of accuracy. 
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number . 



Figure 9 is a plot of J 



and J vs, 
P 



iteration 



Sub-problem number 




Figure 9 



J and J vs. iteration number 
s p 

for missile Intercept problem A. 



After K is increased at the beginning of sub-problem 2, the 
P 

inequality constraint penalty terms become very small. 
This is because the angle of attack is well within its 
admissible region. 

X \ 
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It is evident from Figures 7, 8, and 9 that two 
sub-problems are sufficient to obtain a reasonable solution. 
The overal stopping criterion 

|(Jg*+Jp*)^ - (5.60) 

where i is the sub-problem number is satisfied after the 
third sub-problem. However, at this point the value of Kp 
is 8 which is not large enough to provide desirable 
penalty functionals for the Inequality constraints. It is 
necessary to provide as large a value of Kp as the computer 
will allow for the final sub-problem to insure that the 
minimization is not influenced by the inequality constraint 
penalty terms when the constrained states and controls are 
within their admissible regions. Accordingly, a final sub- 
problem is .performed with Kp = 32. The algorithm is able 
to handle the increase in Kp from 8 to 32 in one step only 
because no constraint boundaries are active in the solution. 
Figure 10 is a plot of the angle of attack expansion 
computed at the end of the last sub-problem. The region of 
admissible angles of attack is shown. 

Figures 11 and 12 are plots of Mach number and 
flight path angle vs. time. In each plot two curves are 
shovmj one is the expansion for the state as computed at the 
end of the last sub-problemj the other is the state trajec- 
tory obtained by numerically integrating the state equations 
with the optimal control expansion. 
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Angle of attack (degrees) 



80 




Figure 10 

Angle of attack vs, time for 
missile Intercept problem A. 




Time (sec.) 



Figure 11 

Mach number vs. time for 
missile Intercept problem A. 
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Figure 12 

Plight path angle vs. time for 
missile intercept problem A. 
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Figure 13 is a plot of relative range X vs. relative 
cross-range Y, As before two curves are presented: one 

represents the expansions of the states as computed at the 
end of the last sub-problem j the other is the state trajec- 
tory obtained by numerically integrating the state equations 
with the optimal control expansion. 



Relative range (ft. x 10^) 



0 2 4 6 8 10 12 




Figure 13 

Relative range vs. relative cross-range 
for missile intercept problem A. 
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Figure is a plot of range X’ vs. cross-range Y* 
where both quantities are obtained by transforming the 
expansions of the states X and Y obtained at the end of the 
last sub-problem from the (X,Y) coordinate system to the 
inertial coordinate system (X’,Y’) fixed at the missile 
launch point (Figure ^10). 




Figure 1^1 

Range vs. cross-range for 
missile intercept problem A. 
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Figure 15 is a plot of load factor vs. time where 



the load factor is given by 

n = I (0 M) (5.61) 

and the states used in equation ( 5 . 61 ) are the state expan- 
sions obtained at the end of the last sub-problem. 




Figure 15 

Load factor vs. time for 
missile intercept problem A. 
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Although load factor is not a state but a function of states, 
the plot is important because it shows that the load factor 
constraints as given by (5.19) are not active. 

It might be suspected that the optimal trajectory 
is a maximum performance turn limited by a constraint 
boundary (angle of attack or load factor) followed by a 
straight line path to intercept. This is not the case. The 
initial turn rate of the missile is small compared to its 
maximum turn rate capability. This is due to the high 
induced drag associated with high angle of attack turns 
which would reduce the missile’s longitudinal acceleration 
capability. Also there is no straight line segment to the 
trajectory although the turn rate of the missile is very 
small at Intercept. 

The first sub-problem was also solved with an 
initial guess for the c vector of: 

all expansion coefficients = 0 
M(T) = 2.5 

e(T) = 0° (5.62 ) 

a(0) = 10° 
a(T) = 0° 

At = 8/20 sec. (T = 8 sec.) 

The sub-problem reached a minimum of 7.5^^ which compared 
favorably with the minimum reached by the first initial 
guess given by equations (5.58). This gives an indication 
that the minimum attained is the global minimum. 
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2. Problem B - 12 Coefficients for each Expansion 
Pour sub-problems were required to solve this 
problem also. Tables 5 and 6 summarize the performance of 
the algorithm. 



sub- 

problem 


e 


r 


K 

P 


I* 


optimization 

strategy** 


C.P.U. 

time 


1 


1.0 X 10“^ 


100 


4 


8 


MMMMMMMM 


31 3411 


2 


0.6? X 10“^ 


100 


6 


2 


PP 


2 '20'' 


3 


O 

X 

O 

1 


100 


8 


2 


PP 


2 '19" 


4 


0.5 X 10“^ 


100 


32 


2 


PP 


0 

C\J 

CVJ 



* number of Iterations required 
** M - MNR method 
P - PNR method 



Table 5 

Weighting factors, optimization strategy, and 
C.P,U. time for missile Intercept problem B. 



sub- 

problem 




J* 


J * 
s 


J * 
P 


1 


11.51 


7.537 


0.2125 X 10“^ 


0.1849 X 10“^ 


2 


9.69 


7.612 


0.1385 X 10“^ 


0.3412 X 10“® 


3 


10.65 


7.628 


0.1513 X 10“^ 


0.3474 X 10"^^ 


4 


10. 31 


7.631 


0.1339 X 10“^ 


0.6139 X 10"^^ 



Table 6 

Components of the minimum augmented measure for 
missile intercept problem B. 
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Figure l6 is a plot of the augmented performance 
measure vs. iteration number and corresponds to Figure 7 
for problem A. 




Figure l6 

Augmented performance measure vs. iteration number 
for missile Intercept problem B. 
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A second failure mode of the MNR method close to the minimum 



is shown. The MNR method produces a divergent in the 
second sub-problem. 

Figure 17 is a plot of range X' vs. cross-range Y' 
obtained in the same manner as in problem A (Figure l4). 




Figure 17 

Range vs. cross-range for 
missile intercept problem B. 
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A comparison of the tables and figures for problems 
A and B show that the optimal control and trajectory have 
not been markedly affected by the increase in the number of 
coefficients from 8 to 12 for each expansion. It is prohib- 
itive in terms of computation time and storage requirements 
to Increase the number of coefficients further. 
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VI. A CLIMB PE RFO PRANCE PROBLEM 

In this section a climb performance problem is solved, 

A supersonic fighter aircraft is to climb from sea level to 
high altitude in minimum time. 

Flight test experience has shown that to climb to 
altitudes above the tropopause in minimum time a supersonic 
fighter must execute a maneuver which typically includes: 

a. a subsonic climb to an altitude near the 
tropopause; 

b. a level or near level acceleration to some 
supersonic Mach number; 

c. a "zoom” climb to the desired altitude trading 
kinetic energy for potential energy. This technique has 
been used extensively in the past decade for establishing 
climb records and in fighter-interceptor tactics to attain 
altitudes higher than the aircraft’s service ceiling for 
short periods of time. 

Several factors contribute to the optimality of this 
type of maneuver. They are: 

a. a fighter's maximum Mach number or "placard" 
limit which arises from dynamic pressure and/or thermal 
limitations and is a function of altitude; 

b. a fighter's transonic drag characteristics; 

c. air temperature variation with altitude; 

d. air density variation with altitude; 

e. turbojet engine maximum thrust variation with 
altitude . 
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Optimization techniques were first applied successfully 
to this problem by Bryson [Refs. 1 and 2], The method of 
steepest descent was used successfully to predict the type 
maneuver described above for a typical supersonic aircraft. 

The epsilon method is applied to the problem herein to 
demonstrate the method's power. A direct comparison of 
methods is not made as the mathematical model used here has 
been Improved considerably over that used in Reference 2. 



A. pr 6 blem formulation 



1. State Equations 

The^ state equations for this problem derived in 
Appendix A are 



M = 



Y = 



gTh g 4 

V " f " -m— 



aM^C, 



gTh sing ^ ^^o^^ 



M 



OTJ 



e 

aMC, - £2^ 



( 6 . 1 ) 

( 6 . 2 ) 



A - aM 
h = 5 — sxn 



( 6 . 3 ) 



The states are Mach number M, vertical flight path angle y» 
and normalized altitude h. The control is angle of attack a. 
Parameters considered constant are the gravitational constant 
g, sea level standard density p^, aircraft wing area S, 
aircraft weight W^, and the altitude of the tropopause under 
standard atmospheric conditions Hj^. These constants are 
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(S.H) 



g = 32.1725 ft. /sec. ^ 

= 0.002378 slugs/ft. ^ 
S = ilOO ft.^ 



= 39,000 lbs. 

Hj^ = 36,089 ft. 



The remaining parameters in equations ( 6 . 1 ) thru 



(6.3) vary with flight and atmospheric conditions. These 
variations are represented in either tabular form or by 
empirical relations. These tabular and empirical relations 
are critical to the problem as they represent the mathemat- 
ical equivalents of the factors a through e listed in the 
introduction to this section. 



With the definition ( 6 . 5 ) Incorporated the state equations 
are 



It is convenient to define the constant 




( 6 . 5 ) 



M cosa - — siny - caoM^Cp. 

a a ' D 



( 6 . 6 ) 



* _ gTh sing 
^ a M 




( 6 . 7 ) 




( 6 . 8 ) 



9 ^ 
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2 , Empirical Relations 

Empirical relations are used for air density ratio a, 
maximum Mach number Mj^, and the speed of sound a as functions 
of normalized altitude h. Air density ratio and normalized 
altitude are defined by 

(6.9) 

( 6 . 10 ) 

These empirical relations which are discussed In Appendix D 
are repeated here for convenience. They are 



-c,h -c^h 

a = e + c^he (6.11) 

= 2.1 - l.le"^*^^ (6'.12) 

a = a (1 - c„h) , h < 1 (6.13) 

o l ’ 

= 971 ft. /sec. , h>l . (6.1^0 

where 

c^ = 1.5^100 (6.15) 

C2 = 1.80ilil5 (6.16) 

c^ = O.illSO (6.17) 

= 0.1331 . (6.18) 



and 



= ^ 



h = JL 



\ 
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3 . Tabular Functions 



In situations where parameter variations cannot be 
adequately represented by empirical formulas, a table of 
values is used. Tables are used for lift and drag coeffi- 
cients as functions of Mach number and angle of attack for 

« 

a typical supersonic fighter. Excerpts from these tables 

are presented in Appendix B. 

The thrust Th appearing in equations (6.6) and (6.7) 

is normalized maxlmiim thrust as it is assumed that since 

the aircraft must climb to altitude in minimum time, its 

power plant will always be operated at fiiaximum thrust. 

Maximum thrust is normalized with respect to sea level 

static maximum thrust Thj^ and is given by 

o 

Th = ^ (6.19) 

o 

where 

ThM = 3^,000 lbs. (6.20) 

Maximum thrust is given as a tabular function of Mach number 
and altitude for the fighter under consideration. 

^ . Performance Measure 

The performance measure for this problem is 



J 




dt 



( 6 . 21 ) 
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5 . Inequality Constraints 

The following state and control Inequality constraints 
are Imposed: 



The maximum Mach number constraint represents the 
"placard" limit. Is a function of altitude and Is given 

by an empirical relation as discussed In Section VI. 2, 



angle of attack slightly above that for maximum lift coeffi- 
cient. The minimum angle of attack Is set slightly below 
that for minimum lift coefficient thus simulating aerodynamic 
stall. 

6. End Conditions 

The Initial conditions are 



0 < M < M, 






( 6 . 22 ) 



- 6 ° = £ ot ^ = 24° 



(6.23) 



The angle of attack constraint Is set at an 



M(0) = M = 0.6 



o 



(6.24) 



y(o) = 0 



(6.25) 



h(0) = h = 0 
o 



(6.26) 



The final condition Is 




(6.27) 
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B. THE EPSILON METHOD FORMULATION 



1 . The Augmented Performance Measure 

Using the penalty functional described in Section III 
for inequality constraints, the augmented performance measure 
is 




+ 



+ 



+ 



+ 




cosa + — siny + caaM^CT^l dt 
a a D J 

m 31- . oaaMC, 1 at 




( 6 . 28 ) 



where e and r are weighting factors, 6 is a constant used to 
make minor adjustments to the admissible regions, and d is 
the midpoint of the admissible angle of attack region; that 
is 



d 



M L 



( 6 . 29 ) 
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The required elements of w are 



«k W ^ ^ " “k'’k\"'^ij [ ff • 



k=l,2, . . . ,K 



w 



K+k 



w 



2K+k 



I T^k “ VlA\ J I e J * 

k=l,2, . . . ,K 



= - 1 



w 



IK 



3K+k " [ ^ 



(rAt)^ , k=l,2,...,K 



w 



4K+k 



r«k"^1^D ^ 

|_a„ - dj (rAt) ® , k-l,2,..,,K 



w^K+i = [(K-l)At]^ 



( 6 . 30 ) 

( 6 . 31 ) 

( 6 . 32 ) 

( 6 . 33 ) 

( 6 . 34 ) 

( 6 . 35 ) 



2 . Functional Expansions, Unknowns, and Partial 
Derivatives 

The state and control expansions are of the same 
form as the problem In Section V and are not shown. The 
elements of the c vector are 



(s^l > ^2 > • • • *^M*^1 *^2* * * * *^M**^1**^2* * * * * *^M* 

( 6 . 36 ) 

d^ ) ^2 » • • • » ^K* ^1 * ^ 
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where the represent the Mach number expansion coeffi- 

cients, the bjjj's represent the vertical flight path angle 
coefficients, the represent the altitude coefficients, 

and the represent the angle of attack coefficients. 

The first and second derivatives of the empirical 
relations (6.11) thru (6.14) are required. These expressions 
are 



d.0 


(6.37) 


— ^ e - [{ 02 - 020 jh)Cj + 


(6.38) 




(6.39) 


= -6.336e-2-‘'h 

dh'^ 


(6.40) 


if " ”^o®7 * h < 1 . . 


(6.41) 


= 0 , h ^-1 


(6.42) 


^ = 0 
dh^ ° * 


(6.43) 


The first and second partlals of the tabular 


functions 



for lift coefficient, drag coefficient, and maximum thrust 
with respect to their Independent variables and the elements 



1 







of w given by equations (6.30) thru (6.35) with respect to 
c are obtained in the same manner as in the problem in 
Section V, 



C. RESULTS 



Two problems were solved. In problem A the aircraft was 
to climb from sea level to 60,000 feet in minimum time. 
Problem B encompassed a s-erles of problems. The results of 
problem A were used as a first guess for the solution of a 
mlnlmum-tlme-to-cllmb profile from sea level to 61,000 feet, 
which in turn was used as a first guess for a climb to 
62,000 feet, etc. In this manner optimal control and state 
trajectories were obtained for minimum-time- to- climb profiles 
from sea level to altitudes from 60,000 to 70,000 feet in 
thousand- foot increments. In both problems 8 coefficients 
for each expansion and ^11 time points (40 time intervals) 
were used. 

1 . Problem A - A Climb from 0 to 60,000 Feet 
The initial guess for the c vector was: 



• all expansion coefficients = 0 
M(T) = 0.9 
y(T) = 45° 
a(0) = 1° 
a(T) = 1° 

At = 120/40 sec. (T = 2 min.) 



(6.44) 



Four sub-problems were required to solve the problem. 
Tables 7 and 8 summarize the performance of the algorithm. 



♦ 



4 i 






I 



sub- 

problem 


e 


r 






optimization 

strategy** 


C.P.U. 


1 


0.2 X lO-^ 


100 


k 


9 


FMMMFFFFF 


3' 33" 


2 


0.1 X 10"^ 


100 


6 


9 


FMMMMMMFF 


3'2iJ" 


3 


0.67 X lO"^ 


100 


8 


8 


MPPPPPPP 


3' 38" 


H 


0.5 X 10"^ 


100 


32 


3 


FFF 


2' 59" 



* number of iterations required 
** M - MNR method 
F - FNR method 



Table 7 

Weighting factors, optimization strategy, and C.P.U. 
time for climb performance problem A. 



sub- 

problem 


Ja* 


J* 






1 


422.5 


129.1 


0.4496 X lO"^ 


0.6855 


2 


442.2 


199.3 


0.2154 X 10"^ 


0.2742 


3 


489. 3 


238.6 


0.1544 X 10"^ 


0.1904 




451.5 


258.0 


0.9677 X 10“^ 


0.8838 X 10"^^ 



Table 8 

Components of the minimum augmented performance measure 
for climb performance problem A. 
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Neither the angle of attack or maximum Mach number 
constraints are active in this problem. The largest Mach 
number attained in the climb is 1.53 at an altitude of 
23,000 feet. Consulting Appendix D, the maximum Mach number 
at this altitude is 1 . 88 , However, since both constrained 
parameters approach their boundaries, a large value of K 

P 

is required to obtain small contributions to the aug- 
mented performance measure. 

Figure 18 is a plot of the augmented performance 
measure vs, iteration number. 




Iteration nurber 
Figure I8 

Augmented performance measure vs. iteration 
number for climb performance problem A. 
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Iterations performed by the FNR method are Indicated by a 
solid line. Iterations performed by the MNR method are 
indicated by a broken line. As observed in the previous 
problem, the FNR method results in excellent terminal 
convergence . The MNR method performs well when the c 
vector is far from optimum as indicated by relatively large 
augmented performance measures. At the commencement of 
sub-problem 1 the MNR method falls presumably because the 
Initial guess for the c vector is too close to optimum. 

The FNR method does not reduce the augmented performance 
measure but manages to salvage the first iteration. On 
Iteration number 2 the opposite occurs. The c vector is 
too far from optimum for the FNR method to work. The MNR 
method comes to the rescue. The same thing occurs at the 
beginning of sub-problem 2, In these two sub-problems the 
use of both methods in combination has allowed the algorithm 
to proceed where the exclusive use of either method by 
itself produces a divergent condition from which the 
algorithm cannot recover. 

Figure 19 is a plot of the performance (final time) 
vs. iteration number. 

Figures 20, 21, and 22 are plots of the states vs. 
time. In each plot two curves are shown: one is the 
expansion for the state as computed at the end of the last 
sub-problem; the second is the state trajectory obtained by 
integrating the state equations with the optimal control 
expansion. 



SUB-PROBLEM 

1 2 3 4 




Figure 19 

Performance measure vs. iteration number 
for climb performance problem A. 
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Mach number 




Time (min.) 



Figure 20 

Mach number vs. time for 
climb performance problem A. 
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Altitude vs. time for 
climb performance problem A. 




Vertical flight path angle vs. 
time for climb performance problem A. 
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Close observation of Figures 20, 21, and 22 reveals a 
trajectory very similar to that described in the introduction 
to this problem. The trajectory begins with a sub-sonic 
climb to an altitude of 33 i 000 feet. The climb angle during 
this portion of the climb reaches a maximum of 27 degrees. 

At this point an acceleration is performed to a Mach number 
of 1.53 with the aircraft in a slight descent. A "zoom" 
climb is then performed to the desired altitude of 60,000 
feet . 

Figure 23 is a plot of the angle of attack expansion 
computed at the end of the last sub-problem. 




“L 

Figure 23 



Angle of attack vs. time 
for climb performance problem A. 
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Initially, the angle of attack is decreasing corresponding 
to the initial acceleration of the aircraft from a starting 
Mach number of 0.6, As the aircraft rotates to climb 
attitude, the angle of attack increases. As the aircraft 
levels off for the supersonic acceleration, the angle of 
attack decreases correspondingly. The angle of attack 
begins to increase again as the "zoom" climb attitude is 
established. A further increase is evident as the aircraft 
slows down in the climb until as the final altitude is 
approached, the aircraft is near stall. The climb was 
completed in 4 minutes and l8v seconds. 

Figure 24 is a plot of altitude vs. range where both 
quantities are obtained by integration. The scales are not 
the same , 




Figure 24 

Altitude vs. range for 
climb performance problem A. 
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The problem was also solved with the maximum Mach 
number restricted to 1.0 throughout the flight regime to 
obtain a comparison of the "zoom” climb technique to a 
totally subsonic climb to 60,000 feet. The aircraft was not 
able to complete the climb. The altitude of 60,000 feet is 
apparently above the service celling of the model. After 
4 minutes and l8 seconds, which was the time required to 
complete the climb by the "zoom” technique, the aircraft 
was passing ^ 13,000 feet and climbing very slowly. 

2 . Problem B - Optimum Climbs to Altitudes 
from 60,000 to 70 » 000 Feet 

Table 9 depicts the results for each sub-problem as 
minimum time-to-climb profiles are generated for final 
altitudes of 60,000 to 70,000 feet in thousand-foot increments. 
For each sub-problem the results of the previous sub-problem 
were used as a first guess for the new trajectory. The 
stable behavior of the FNR method near the minimum is respon- 
sible for the ability of the algorithm to generate neighboring 
optimal trajectories. Thus, in effect, ten problems were 
solved with minimum effort by taking advantage of the results 
of each problem in turn. If, however, the change in end 
conditions is too large, the MNR method may be required to 
start the new sub-problem. This is the case in sub-problems 
11 thru 1^. 

Figure 25 is a plot of altitude vs. range where 
both quantities are obtained by numerical integration showing 
the climb trajectories for several of the final altitudes. 
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sub- 

problem 


nunher of 

iterations 

required 


optimization 
strategy * 


final 
altitude 
( feet ) 


time to 

climb 

(J*) 


C.P.U. 

time 


l\ 


3 


PPP 


60,000 


4 '18" 


2 '59" 
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PPPP 


61,000 


4 '24" 


3'11" 


6 


5 


PPPFF 


62,000 


4' 31" 


3' 16 " 


7 


12 


PPPFFPPFPFPP 


63,000 


4' 36 " 


5'02" 


8 




PPPP 


64,000 


4' 4l" 


3 '05" 


9 


3 


MPP 


65,000 


4'44" 


2 '52" 


10 


5 


PPFPP 


66,000 


4 '54" 


3'29" 


11 




MPPP 


67,000 


4' 57" 


3'04" 


12 


3 


MPP 


68,000 


5 '00" 


2'50" 


13 


4 


MPPP 


69,000 


5 ’03" 


2 ’55" 


lH 


4 


MPPP 


70,000 ; 


5'09" 


2 '53" 



* p _ FNR method 
M - MNR method 



Table 9 

Results of climb performance problem B. 
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Figure 25 

Altitude vs. range for 
climb performance problem B, 



The climb profiles shown in Figure 25 reveal the following 
characteristics: 

a. The sub-sonic climb profiles are identical for 
all final altitudes for the initial portion of the climb; 

b, as the final altitude increases, the altitude at 
which the aircraft levels off for the supersonic acceleration 
increases by approximately 15 to l8 percent of the final 
altitude ; 
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c. the aircraft performcs a diving maneuver to 
transit the transonic region with the maximum dive angle 
(13°) reached In the climb to 60,000 feet; 

d. as the final altitude Increases, the aircraft 
performs the supersonic acceleration at higher altitudes 
with less altitude loss In acceleration. The results are 
not In agreement with standard practice In which the 
accelerations are generally performed In level flight at 
the tropopause. The results are in agreement with the 
results of Bryson [Ref. 2] which clearly show the dive 
associated with transiting the transonic region. Bryson’s 
results were obtained by the method of steepest descent. 
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VII. AN AIR COMBAT MANEUVERING PROBLEM 



In this section the turning performance of a supersonic 
fighter Is considered. First, the basic aircraft limitations 
on maneuvering flight are reviewed. Second, turning perfor- 
mance In the horizontal plane Is reviewed from a theoretical 

point of view. The "corner” velocity concept familiar to 

% 

fighter pilots Is presented. Third, turning performance In 
three-dimensions Is discussed. Finally, a three-dimensional 
problem Is solved In which a supersonic fighter Is required 
to execute a l80® course reversal In minimum time with the 
Initial and final altitudes specified. 

A. THEORETICAL TURNING PERFORMANCE 

1 . Aircraft Performance and Maneuvering Limitations 

A tactical fighter must be highly maneuverable. An 
Important consideration In maneuverability Is the ability 
of the fighter to turn. Two basic performance criteria In 
turning performance are : 

a. radius of turn, and 

b. rate of turn. 

In air combat situations It Is often desirable to perform 
a turn so that the aircraft's radius of turn (curvature) Is 
minimized or the aircraft's rate of turn Is maximized. The 
ability of a fighter to minimize radius of turn or maximize 
rate of turn Is limited by: 



a. maximum thrust, 

b. aerodynamic stall, or 

c. maximum allowable load factor. 

2. Turns in the Horizontal Plane 

If an aircraft is restricted to move in a horizontal 
plane only, turning performance is easily analyzed. Using 
the assumptions given in Appendix A and the added assumption 
that 



T sin a << L, 



(7.1) 



equations for lift coefficient, radius of turn, rate of 
turn, and the thrust required to maintain level flight at 
constant velocity are easily derived and well known. They 
are 






•2W^n 

pSv^ 



R = 



g(n^-l) 



(7.2) 



(7.3) 






g(n^-l) 

V 




( 7 .^ 4 ) 



and 

I 

T = v^ Cj. . (7.5) ■ 

2 cos a 
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where 



T = thrust, 
a = angle of attack, 

L = lift, 

Cj^ = lift coefficient, 

W = aircraft weight , 
n = load factor 
p = air density, 

S = wing area, 

V = aircraft velocity, 
g = gravitational constant, 

Cp = drag coefficient, 

R = radius of turn, and 

^ - horizontal flight path angle (heading angle) . 



Figure 26 is a plot of lines of constant thrust, constant 
radius of turn, constant rate of turn, and maximum lift 
coefficient on velocity-load factor CV-n) diagrams. 



constant 
rate of turn r 



constant 

radius of turn R 



constant 
thrust T 




Figure 26 

Velocity-Load Factor Diagrams 
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In Figure 26(c) the constant thrust lines indicate the 

thrust required to maintain a steady state turn at a specific 

load factor and velocity. The corner velocity v is defined 

c 

as that velocity at which an aircraft is capable of operating 
at maximum lift coefficient C- and maximum structural load 
factor nj^ at the same time. This is the velocity which 
produces minimum radius of turn and maximum rate of turn as 
can be seen from Figures 26(a) and 26(b). The corner 
velocity can only be maintained in the steady state if the 
aircraft has enough thrust available to allow the maximum 
thrust curve to pass above the corner created by the 
Ct - n.« boundary intersection. If sufficient thrust is 
not available to allow the corner velocity to be maintained 
in the steady state, which is typically the case, the air- 
craft must either degrade its turning performance by moving 
off the boundary intersection until the maximum thrust curve 
is encountered, or sacrifice altitude. In this case the 
velocity for maximum rate of turn is larger than the velocity 
for minimum radius of. turn. 

As can be seen from the previous discussion, 
optimization techniques are not required to analyze turns in the 
horizontal plane. This arena, however, is an excellent 
place to test an optimization technique which is being 
considered for use in solving more complicated problems 
Involving three-dimensional maneuvering. With this in mind, 
two problems involving turns in the horizontal plane were 
solved by the epsilon method and the answers compared to 
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theoretical results. In one problem an aircraft v;as 
required to perform a horizontal turn with minimum radius. 

In a second problem the aircraft was required to turn 
through a given heading change in minimum time which is 
equivalent to maximizing rate of turn. The aircraft was 
given a large thrust capability so that the turns were not 
thrust limited. The mathematical model used in the problems 
is given in Appendix A. The model is an accurate three- 
degree-of- freedom model of the aircraft’s motion with 
maneuvering limitations included. Prom the previous 
discussion the aircraft should have performed the turn 
in both cases at the corner velocity where 



V. 



. r 






[pS(Cj^ tana^ + ) 



(7.6) 



"M 



M 



In this equation a is the angle of attack for maximum lift 

s 

coefficient and C„ is the drag coefficient for maximum lift 
coefficient. The epsilon method solved both problems 
successfully. In each case the optimum trajectory Involved 
a turn at the corner velocity. Thus, the ability of the 
second-order epsilon method to handle problems of this type 
was demonstrated. 

3. Turns in Three Dimensions 

The analysis of turning performance in three 
dimensions is quite complicated. In this regime modern 
optimization techniques are the only method of solving 
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meaningful problems. Even optimization techniques are apt 
to have a difficult time with three-dimensional maneuvering 
problems using realistic models of the aircraft's motion 
because of large computer time and storage requirements. In 
this section the epsilon method is used to solve an impor- 
tant three-dimensional maneuvering problem often encountered 
in air combat . 

In many air combat situations a fighter pilot is 
faced with the requirement to reverse his course as rapidly 
as possible. Generally, the pilot has in mind a specific 
altitude at which he would like to complete the maneuver 
which is dictated by his desire to track an enemy aircraft 
or perform some attacking maneuver. With this in mind, the 
problem posed for solution by the epsilon method involves a 
supersonic fighter which is required to execute a l80® 
reversal in minimum time. The aircraft must begin the 
maneuver in level flight and recover in level flight at the 
entry altitude. The accepted maneuvers used to accomplish 
this task developed over years of combat experience are the 
high-speed yo-yo and the low-speed yo-yo maneuvers. If the 
aircraft begins a reversal at a flight speed higher than 
its corner velocity a high-speed yo-yo is called for and 
vice-versa. A high speed yo-yo consists of a climbing turn 
followed by a descending turn. A low speed yo-yo consists 
of a descending turn followed by a climbing turn. If the 
aircraft begins a reversal at its corner velocity, a level 
turn is called for. The purpose of applying the epsilon 
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method to this problem is to either confirm or challenge 
the effectiveness of these experimentally developed 
maneuvers by the use of an optimization technique. The 
assumptions applied to the problem, the coordinate system 
and nomenclature used, and the derivation of the equations 
of motion are presented in Appendix A. 



B. PROBLEM FORMULATION 
1 . State Equations 

The state equations derived in Appendix A are 



M 



Thcosa - I slny 



gpS,a 

2W ‘‘ 
e 



( 7 . 7 ) 



• £ nsln4>. 

^ a Mcosy 



( 7 . 8 ) 



• £ ncos4) _ £ cosy gx 

h = ^ siny . (7.10) 

“l 

The states are Mach number M, horizontal flight path angle 
\lif vertical flight path angle y, and normalized altitude h. 
The controls are bank angle (j), normalized thrust Th, angle 
of attack a, and load factor n. 

In addition to the four state equations, an equality 
constraint which must be satisfied is 



0 



> Thslnc. . # 
e e 



M Cj^ - n 



( 7 . 11 ) 



120 



I 



This equation is a result of the definition of load factor 
n and is derived in Appendix A. It is possible to use 
equation (7.11) to eliminate the load factor control from 
the state equations, but this Is not desirable for two 
reasons: first, the resulting state equations would be 

further complicated thus Increasing the analytic v;orkload 
required to compute first and second partial derivatives; 
second, the incorporation of the Important load factor 
Inequality constraint would be unnecessarily complicated. 

Parameters considered constant for this problem are 
the gravitational constant g, aircraft maximum thrust Thj^, 
aircraft weight W^, speed of sound a, base altitude Hj^, air 
density p, and aircraft wing area S. These constants are 



g = 32,1725 ft. /sec 
Thj^^ = 21,000 lbs. 

Wg = 39,000 lbs. 
a = 1077.8 ft. /sec. 



2 



( 7 . 12 ) 



= 10,000 ft. 

p = 0.001756 slugs/ft. ^ 

S = 400 ft.^ 



It Is convenient to define the constants 




(7.13) 




(7.14) 



°2 “ W a 
e 




(7.15) 
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( 7 . 16 ) 



e 







(7.17) 



With these simplifying definitions incorporated the state 
equations are 



2. Tabular Functions 

Tables are used for lift and drag coefficient as 
functions of Mach number and angle of attack for a typical 
supersonic fighter. Excerpts from these tables are provided 
in Appendix B. 




( 7 . 18 ) 




(7.19) 




( 7 . 20 ) 




( 7 . 21 ) 



and the additional equality constraint is 




( 7 . 22 ) 
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Performance Measure 



The performance measure for this problem is 







J = / dt . 


(7.23) 


Inequality Constraints 




The controls must satisfy 




0 < Th < 1 


(7.24) 


0 ^ a < ttjyj = 24° 


(7.25) 


0^n <_nj^= 6.5 g’s 


(7.26) 


0 <_ (j> <_ IT . 


(7.27) 



The minimum allowable normalized thrust, angle of attack, 
and load factor are approximated by zero as these constraints 
are not anticipated to be active, A zero value of the 
lower bound simplifies the associated penalty term in the 
augmented performance measure. The maximum angle of attack 
is fixed at a value slightly higher than the angle of attack 
for maximum lift coefficient as given in the tabular data 
thus simulating aerodynamic stall. The structural load 
factor upper bound is 6.5 g's, a standard value from fighter 
aircraft operational limitations. The bank angle constraints 



123 



1 



were required to keep the algorithm from generating bank 
angles greater than l80®. 



5. End Conditions 

The initial conditions are 

M(0) = 0.9 

ip(0) = 0 

y(o) = 0 

h(0) = h = 15,000 ft. 

o * 

The final conditions are 
. \|>(T) = TT 
y(t) = 0 
h(T) = h^ . 



( 7 . 28 ) 

( 7 . 29 ) 

( 7 . 30 ) 

( 7 . 31 ) 

( 7 . 32 ) 

( 7 . 33 ) 
( 7 . 3 ^) 



C. THE EPSILON METHOD FORMULATION 

1. The Augmented Performance Measure 

Using the penalty functionals described in Section 

s 

III for the inequality constraints, the augmented perform- 
ance measure is 



i \ 



12k 



dt 



-/ 



+ I^M - C 2 Thcosa + c^siny + dt 



^ rfj fi - « 

* r/ fv - =1 - =1 « 



. 1 f [: an • „12 

* rj[, h ■ H^; 



c^jThsina + c^M^Cj^ - nj 



dt 



(7.35) 






where e and r are weighting factors, and 6 is a constant 
used to make minor adjustments to the admissible regions. 

The required elements of v/ are 

'"k " ■ C2l\cosotj^ + c^slnyj^ + , k^l,2, . . . ,K (7.36) 



_ r *. TAtlJ^ , ^ 

'^K+k L'^'k ■ °1 Mj^cosYj^J L ej > k 1,2,...,K 



(7.37) 
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'^SK+k 


[\ - f [f]^ 


, k=l,2,...,K 


(7. 39) 


'^^K+k 


[c^Th^sino^ + • 




(7.40) 


^5K+k " 1 


-2Th, IK , 

- ij P(rit)’« 
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(7.41) 


'^bK+k " 1 


e - ■ 


k=l,2, . . . ,K 


(7.42) 


'^TK+k " [ 




k=l,2,...,K 


(7.43) 


'^SK+k " 1 




k=l,2, . . . ,K 


(7.44) 



W 9 K +1 = C(K-l)At]^ , (7.45) 

2, Functional Expansions, Unknowns, and 
Partial Derivatives 

The state and control expansions are of the same 
form as In previous problems and are not shown. The elements 
of the c vector are 



T 

c — ( 3.^ j ^ 2 9 • • • 9 ** ^ ^ ^2 ^ * * * ^ ^ 

** * * * ^2 ^ * * * * ^ ®2 ***** *^1*^2^ * * * * ^ 

(7.^6) 
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where the 3 -jjj's represent the Mach number expansion coeffi- 
cients, the bj^’s represent the horizontal flight path 
coefficients, the represent the vertical flight path 

angle coefficients, the d *s represent the altitude coeffi- 
dents, the represent the bank angle coefficients, the 

f^'s represent the thrust coefficients, the represent 

the angle of attack coefficients, and the h^*s represent 
the load factor coefficients. 

The first and second partial derivatives of the 
tabular functions for lift coefficient and drag coefficient 
with respect to their Independent variables and equations 
(7.36) thru (7.^5) with respect to the c vector are obtained 
in the same manner as in previous problems. 

D. RESULTS 

Three problems were solved. In problem A the aircraft 
must perform the I80® reversal in minimum time starting 
from an initial Mach number of 0 . 9 . In problem B the air- 
craft must perform the reversal starting from its corner 
Mach number which from equation (7.6) is O.708. In problem 
C the aircraft must perform the reversal starting from an 
initial Mach number of 0 . 5 . In all problems 8 coefficients 
for each expansion and 21 time points (20 time Intervals) 
are used. 

1 . Problem A 

Since the initial Mach number is above the corner 
Mach number, a high-speed yo-yo maneuver is called for by 
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accepted tactics. With this in mind an initial guess of 
the c vector was made to reflect this type of maneuver. 
Accordingly, the following coefficients were given non-zero 
initial values: 

a^ = -0.309 
c^ = 0.52^1 

d, = 0.333 

^ (7.^) 

e^ = l.Oil? 

g^ = 0.262 

h^ = 3.000 



The remaining initial values for the c vector were: 



remaining expansion coefficients = 0 
M(T) = 0.9 
<j)(0) = 0° 

<j)(T) = 0° 

Th(0) = 0.88 (30,000 lbs.) 

Th(T) = 0.88 (30,000 lbs.) 
a(0) = 5° 
a(T) = 5 ° 
n(0) = 1.0 
n(T) = 1.0 

At = 12/20 sec. (T = 12 sec.) 



( 7 .^ 18 ) 



The problem was solved in six sub-problems. It took 
17.65 seconds to complete the turn. Figures 27 thru 30 
are plots of the control expansions as computed at the end 
of the last sub-problem. 

Prom Figure 28 it is seen that the maximum thrust 
constraint is active for the first 10 seconds of the turn. 
At t « 6 seconds there is a short period in which the 
thrust is slightly inadmissible. This is due to the use 
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Thrust (lbs. x 10'^) , Bank angle (degrees) 




Figure 27 

Bank angle vs. time for 
turning performance problem A. 




Figure 28 

Thrust vs. time for 

turning performance problem A. 
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Figure 29 

Angle of attack vs. time 

for turning performance problem A. 




Figure 30 

Load factor vs. time for 
turning performance problem A. 
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of the factor 6 = 1,03 in the inequality constraint penalty 
terms in equation (7.35) which has the effect of slightly 
Increasing the size of the admissible region. This is 
desirable, however, as the epsilon method generates only an 
approximation to the optimal control from which the true 
optimal control must be deduced. It is easier to recognize 
an optimal control expansion which is on a constraint 
boundary with the factor 6 included. As shown in Figure 6 
in Section V, 6 = 1.03 is the proper choice for a final 
Kp = 30. During the last portion of the turn, the aircraft 
is operated at the angle of attack for maximum lift coeffi- 
cient (approximately 22°), The bank angle and load factor 
constraints are not active during the maneuver. 

Figures 31 thru 3^ are plots of the states vs. time. 
In each plot two curves are shown: one is the expansion for 

the states as computed at the end of the last sub-problem; 
the second is the state trajectory obtained by numerically 
integrating the state equations with the optimal control 
expansions. An observation of these plots reveals that a 
high-speed yo-yo maneuver is performed although the altitude 
excursions are not, as large as this author expected. The 
optimization procedure settles on a nearly level hard turn 
at high load factors, steep bank angles, and maximum thrust 
for the majority of the turn. When the maximum thrust 
boundary is not active, the aircraft flies at the angle of 
attack for maximum lift coefficient. 
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Figure 31 

Mach number vs. time 

for turning performance problem A. 




Figure 32 

Horizontal flight path angle vs. time 
for turning performance problem A. 




Figure 33 

Vertical flight path angle vs. time 
for turning performance problem A. 
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Figure 3^ 

Altitude vs. time for 
turning performance problem A. 

2 . Problems B and C 

Figures 35, 36, and 37 are plots of cross-range vs. 
range, altitude vs. cross-range, and altitude vs. range 
obtained by integrating the equations of motion with the 
optimal control expansions found in problems A, B, and C. 

An observation of Figures 35, 36, and 37 reveals that the 
expected maneuvers are performed for each initial Mach 
number. In problem B the aircraft performs an essentially 
level turn from an initial Mach number equal to its corner 
Mach number at this altitude. In problem C the aircraft 
performs a low-speed yo-yo maneuver losing a maximum of 
800 feet after 90° of turn from an initial Mach number 
below its corner Mach number. In problems A and C, however, 
the aircraft does not go through as much of an altitude 
excursion as anticipated by the author. Since in fighter 
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Figure 35 

Cross-range vs, range for 
turning performance problems A, 



B, and C. 
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Altitude vs. cross-range for 
turning performance problems A, 




Figure 37 

Altitude vs. range for 
turning performance problem A, B, 



, and C , 



and C. 
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tactics, however, ther^ are no rules on the amount of 
altitude which should be gained or lost in a yo-yo maneuver 
a quantitative evaluation of the results is purely subjec- 
tive. The important result is that the optimization tech- 
nique did require the aircraft to perform the high-speed 
and low-speed yo-yo maneuvers predicted by accepted tactics 
The accepted tactics are, therefore, qualitatively correct. 

The optimal times required to complete the turn for 
each of the three problems are given in Table 10. 





Optimal time 
for reversal 


problem A 
problem B 
problem C 


17.6 sec . 
20.8 sec. 
2^.9 sec. 



Table 10 



VIII. SUMMARY AND CONCLUSIONS 



A number of realistic problems in aircraft and missile 
performance optimization have been solved by the use of a 
second-order epsilon method. The mathematical models have 
portrayed aircraft and missile dynamics in an accurate 
manner with particular emphasis placed on the modeling of 
performance and maneuvering limitations. 

The state and control inequality constraints generated 
by these limitations have been handled by a new computa- 
tionally superior penalty functional. Three desirable 
theoretical properties of this penalty functional have been 
shown. 

A full Newton-Raphson method for minimizing the aug- 
mented performance measure, has been shown to be computation- 
ally feasible and superior in certain situations to the 
"modified" Newton-Raphson method proposed elsewhere. 

The following observations are significant with respect 
to the second-order epsilon method. 

a. The MNR method is relatively Insensitive to the 
starting values of the unknowns c. The FNR method diverges 
for starting values of c which are far from optimum. 

b. Once c is close to optimum, the FNR method 
converges rapidly whereas typically the MNR method either 
diverges, oscillates, or converges slowly at best. 

c. In relatively simple problems the MNR method is 
capable of obtaining a solution by Itself. In more difficult 
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problems such as those solved in this dissertation, a 
combination of the two methods is required. Typically, the 
most effective procedure involves using the MNR method 
initially followed by the FNR method v/hen successive itera- 
tions yield "small" Improvements in the augmented perform- 
ance measure. In other rare cases where the initial guess 
for the c vector is close to optimum, the FNR method must 
be used initially. 

d. The power of the FNR method close to the minimum 
can be used to advantage to obtain with minimum effort 
optimal control and state trajectories for problems with 
neighboring end conditions by using the solution to a basic 
problem as a first guess for the new problem. 

The solutions to the problems solved have a number of 
applications. In the missile intercept problem (Section V) 
minimum-time optimal trajectories may be used as a basis for 
comparison with the performance of more practical sub- 
optimal controllers such as proportional navigation for a 
short range air-to-air missile . In the three-dimensional 
turn-reversal problem the qualitative optimality of an 
experimentally developed air combat maneuver is shown for 
the first time. A significant lesson to be learned from 
the results of this problem is the' Importance of thrust in 
comparison to lift coefficient in the maneuvering capability 
of a fighter aircraft. Thus, an optimization method of the 
type used in this work applied to realistic performance 
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problems has application In the evaluation of tradeoffs In 
the design of future flight vehicles. 



APPENDIX A 



MATHEMATICAL MODELS 

In this Appendix the mathematical models used In the 
problems are derived. In Section A.l the basic equations 
of motion of an aircraft In three dimensions are derived 
under appropriate assumptions. This model Is used In the 
problem solved In Section VII. In Section A. 2 the aircraft 
Is restricted to move In the horizontal plane only and the 
appropriate adjustments are made to the three-dimensional 
model. In Section A. 3 the aircraft Is restricted to move 
In the vertical plane only and the appropriate adjustments 
are made to the three-dimensional model. This model Is 
used In the problem solved In Section VI. In Section A.^ 
the mathematical model for the missile Intercept problem 
solved In Section V Is derived. 

1 . The Mathematical Model for an Aircraft Maneuvering 
In Three Dimensions 

In this section the basic three-degree-of-freedom 
equations of motion of an aircraft are derived. The 
assumptions are 

a. the earth Is flat, 

b. the aircraft Is a point mass, 

c. the mass of the aircraft Is a constant, 

d. the aircraft Is always In balanced flight, 

e. the aircraft rolls about Its velocity vector, 

f. acceleration due to gravity Is a constant. 



and 
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The coordinate system and notation are presented below. 




Figure 39 

Aircraft Coordinate System 



Three axis systems. are drawn in Figure 39. They are: 



a. (X, Y, H) 

b. (X', y, Z') 



(x, y, z) 



a fixed inertial axis system; 

a non-rotating axis system fixed 
to the center of mass of the 
aircraft; 

a rotating axis system fixed to 
the center of mass of the aircraft; 
the x axis is oriented in the 
direction of the aircraft's velocity 
vector; the y axis points out the 
right wing. 



The attitude of the aircraft is given by four angles as 
follows : 

a. a rotation in the X'Y’ (horizontal) plane; 

b. a rotation y in the xZ’ (vertical) plane; 

c. a rotation <J> in the yz plane; 

d. a rotation a in the xz plane. 

The three angles y, and are the Euler angles (Ref. 19)- 
The angle a is the angle of attack of the aircraft using 
the thrust line as a reference. The remaining notation is 

a. forces; 

L = lift, 

D = drag, 

T = thrust, 

Th = normalized thrust, 

Wg = gross weight, 

b. angles; 

a = angle of attack of the thrust line 
measured in the xz plane, 

a = angle of attack for maximum lift coefficient, 

o 

y = vertical flight path angle measured in the 
xZ* plane, 

(p = bank angle measured in the yz plane, 

ip = horizontal flight path angle measured in 

the X’Y' plane, 

c. rates; 

p = roll rate measured in the yz plane, 
q = pitch rate measured in the xz plane, 
r = yaw rate measured in the xy plane, 

0 ) = angular velocity of the xyz system with 

respect to the X'Y'Z' system, 

d. other parameters; 

V = velocity, 

m = mass, 

g = gravitational constant. 
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p = air density, 

S ®= aircraft wing area, 

= lift coefficient, 

Cp = drag coefficient, 

H = altitude, 

h = normalized altitude, 

R = radius of turn, 

M = Mach number, 

M = Corner Mach number, 

a = speed of sound, 

= corner velocity, 
n = load factor, 

e = unit vector (with appropriate subscript 

indicating direction), 

e. subscripts; 

M = maximum value, 

L = minimum value. 

The equations of motion are derived following the 



methods used in Reference 19. Starting with Newton’s 
second law 



The aircraft velocity and acceleration may be written 




(A.l) 



where ^ is defined as the time derivative in the xyz system 



V = ve 



(A. 2) 




(A. 3) 



where 
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(A.H) 




and 



v®x ^ ^ Y * (A. 5) 

The angular velocity of the xyz system v;ith respect to the 
non-rotating frame X'Y'Z' is given by 



<*) = Pe^ + qCy + re^ . (A. 6) 

At this point, relations between the angular rates p, q, 
and r and the angular rates of change of the Euler angles 
are required. These relations are purely trigonometric in 
nature and are derived in Reference 19. In matrix notation 
they are 



p 




1 0 -slnY 




' 4> ' 


q 


= 


0 coS((> sin4> cosy 




• 

Y 


r 




0 -sin(}) coS(() COSY 







Substituting equations (A. 7) into equation (A. 6), we obtain 

• • • • 

w = ((j) - i|JsinY)e„ + (ycos(}> + \{»sin(}> cosY)e„ 

(A. 8) 

+ (\|)Cos4> COSY - Ysin(}>)e . 
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Using equations (A. 8) and (A. 2) the product 



w X V = (ii)cos({) COSY - Ysin4»)ve, - (li^sincf) cosy + Ycos4>)ve 

(A. 9) 

is formed. Isolating thrust and v/eight components in the 
xyz system. 



T = Tcosa , 

Ty = 0 , 

= -Tsina , 

Wp = -W sinY , 

W = W COSY sin(f) , 

W = WgCOSY cos4) , 

equation (A.l) may be written in component form as 



(A. 10) 
(A. 11) 
(A. 12) 
(A. 13) 

(A. 14) 

(A. 15) 



Tcosa - D - WgSinY = mv , (A. 16) 

• • 

¥ COSY sincf) = mv(i/jcos({) cosy ~ Ysln4)) , (A. 17) 

• ^ 

WgCOSY cos4» - Tsina - L = -mv(i|jsin4» cosy + ycos4)).(A.18) 



The equations (A.l6) thru (A.l8) are the basic equations of 
motion. To apply optimization methods, it is desirable to 



transform these equations into state variable format. 

First, lift and drag coefficients are defined by the 
expressions 

(A. 19) 
(A. 20) 

Second, it is convenient to introduce the normalizing 
expressions 



2 

L = Cj^J^pv S , 
D = Cjj%pv^S . 



V = aM , 



(A. 21) 



T = Thj^Th . 



(A. 22) 



Substituting equations (A. 19) thru (A. 22) into equations 
(A.16) thru (A.18) along with the expression 

Wg = mg , (A. 23) 

we obtain 



o W 

Thj^^Thcosa -3gCpp(aM)^S - W^siny = aM 


(A.2i|) 


W . 

W cosysincf) = — ^ aM(if)eos(J) cosy - ysin4>) 

^ s 


(A. 25) 


WgCosy cos(}) - Thj^Thsina - %Cj^p(aM)^S 


(A. 26) 


W . . 

= aM(iJjsin(j) cosy + ycoscj)). 

o 




1^17 

• 


1 



h 



I 









Solving equations (A. 2^4) thru (A. 26) for M, and y yields 
the state variable format 



M = 


gThj^ 

W a 
e 


Thcosa - ■§ 
a 


slnY - 

e 


(A. 27) 


i = 


gThn 


Thsina sin<{) 


, KpSa 


(A. 28) 


W a 
e 


Mcosy 


2Wg cosy 


• 

Y = 


gThj,^ 

W a 


Thsina cos<() 
M 


gpSa g cosy 

2W a M 

P 


(A. 29) 



In addition the position of the aircraft (center of mass 
location) may be required from some fixed reference point. 
To this end three additional state equations are 



X = aMcosy cosif) , (A. 30) 
Y = -aMcosy sini|^ , (A. 31) 
ft = aMsiny . (A. 32) 



It is convenient, also, to define the load factor n as 



n = 



- TOTAL LIFTING FORCE 



WEIGHT 



L + Tsina 
W. 



^ Thslna + m\ 



(A. 33) 
(A. 3^) 

(A. 35) 



With equations (A. 35) incorporated into equations (A. 27) 
thru (A. 29 ) the state equations are 



M = -rt — Thcosa 
W a 




(A. 36 ) 



: ^ £ nsln4> 



a Mcosy 



(A. 37) 



‘ _ g ncos4> g cosy 
^ a M ~ a M 



(A. 38 ) 



The mathematical model for the three-dimensional 



reversal problem solved in Section VII includes the state 
equations (A. 36 ), (A.37)> (A. 38 ), and (A. 32). In addition, 
equation (A. 35) must be satisfied. This equation is written 
in the form 



The states are Mach number M, horizontal flight path angle 
ipf and vertical flight path angle y. The controls are 
bank angle ({>, normalized thrust Th, angle of attack a, 
and load factor n; The purpose of introducing load factor 
as an independent control through equation (A. 35) vice using 
the state equations (A. 27) thru (A. 29) is to simplify the 
state equations and the incorporation of the structural load 




n . 



(A. 39) 



e 



e 



factor constraint. 



The following inequality constraints are imposed 



on the controls: 



0 < Th < 1 , 



(A.ilO) 



0 < a < a. 



M ’ 



(A. 1^1) 



0 < n < n, 



M 



(A.il2) 



The lift and drag coefficients are given as tabular 



functions of Mach number and angle of attack. Reynold's 
number effects are neglected. The parameters considered 
constant for the problem in Section VII are the gravitational 
constant g, maximum thrust Thj^, aircraft gross weight 
the speed of sound a, air density p, and wind area S. 

2 . The Mathematical Model for an Aircraft Maneuvering 
in the Horizontal Plane 

In this Section the state equations (A. 36 ) thru 
(A. 38 ) are applied to an aircraft restricted to maneuver 
in a horizontal plane. The appropriate assumptions are 



Y = 0 



(A.il3) 



(A.kk) 



H = H 



U.^5) 



o 



\ 
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Applying equations (A. ^3) thru (A. ^5) to equation (A. 38 ), 
we obtain 



n = 



COS<() * 






Substituting equations (A. ^3) and (A. 46) into equations 
(A. 36 ) and (A. 37), we may write the state equations as 



M = 



gTh 



M 



W a 
e 



Thcosa - 



gpSa 

2W “ ’ 

e 



(A. 47) 



i = S^ni . 

^ aM 



(A. 48) 



The mathematical model for the two dimensional 
minimum time and minimum radius of turn problems referred 
to in Section VII includes the state equations (A. 47) and 
(A. 48). In addition, equation (A. 35) must be satisfied. 
Using equation (A. 46), equation (A. 35) is written in the 
form 



0 = 



Th 

“ 



M 



Thsina cos<J> + 



pSa' 

2W 



Cj^M cosij) - 1 



(A. 49 ) 



The states are Mach number M, and horizontal flight path 
angle The controls are bank angle 4>, angle of attack a, 
and thrust Th. It is possible to eliminate one control by 
substituting equation (A. 49) into equation (A. 48). The use 
of equation (A. 49) as an additional equality constraint is 
preferred, however, because the state equations are simpler 
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and the required control inequality constraints are simpler 
to incorporate . 

The following Inequality constraints are imposed on 
the controls: 

0 £ Th < 1 , (A. 50) 

0 1 1 4>iyi = cos"^[^] , (A. 51) 

1 1 • (A. 52) 

The lift and drag coefficients are given as tabular 
functions of Mach number and angle of attack. The parameters 
considered constant for the problem are the same as those 
listed for the three dimensional model described in Section 

A.l. 

3 . The Mathematical Model for an Aircraft 
Maneuvering in the Vertical Plane 

In this section the state equations (A. 27) thru (A. 29) 
are applied to an aircraft restricted to maneuver in the 
vertical plane. The appropriate assumptions are 

<|) = 0 (A. 53) 

and 

i = 0 . (A. 54) 
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Substituting equations (A. 53) and (A. 5^) into equations 
(A. 28 ) and (A. 29), we may write the state equations as 



M = ; Thcosa 
w_a 



- ^ siny 



_ ££Sa j^2 

2W " 



C = Thslna . gpSa £ cosy 

’ W^a M 2W L "a M 

e e 



(A. 55) 
(A. 56) 



The mathematical model for the tv;o dimensional climb 
performance problem solved in Section VI includes the state 
equations (A. 55) and (A. 56 ) along with state equation (A. 32). 
It is convenient, however, to introduce the following 
relations into the state equations: 



Q 

11 


(A. 57) 


Thj/j 

Th - 

W 

e 


(A. 58 ) 


^ - h7 


(A. 59) 



where 



and 



Q = 

Po = 

Th = 

«L = 
h = 



density ratio, 
standard sea level density, 
normalized maximum thrust, 
tropopause altitude, 
normalized altitude. 
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Substituting equations (A. 57) thru (A. 59) into equations 
(A. 55) and (A. 56), we obtain the revised state equations 



M 






R SP^Sa 

cosa - “ sina - — oM C,, 

a 2V/ D 



gP^Sa 



Y = 



gTh sina . o 
a M 2W, 



oMC, - 

L aM 



i aM . 

h = sin Y • 



(A. 60) 
(A. 61) 



(A. 62) 



The states are Mach number M and vertical flight path 
angle y. The control is angle of attack a. 

The following inequality constraints are imposed on 
the states and controls: 



0 < M < M„ 

— — M 



(A. 63) 



a 



min 



< a < 



a 



M 



(A.6i|) 



Thrust (Th) represents normalized maximum thrust for 
the problem in Section VI. This is given as a tabular 
function of Mach number and altitude. The lift and drag 
coefficients are given as tabular functions of Mach number 
and angle of attack. 

Empirical relations are used for density ratio a, 
speed of sound a, and maximum Mach number Mj^ as functions of 
altitude. These relations are given in Appendix D. 



15 ^ 





















1 




The parameters considered constant for the problem 
are the gravitational constant g, standard sea level density 
p^, wing area S, gross v/eight V/^, and tropopause altitude 

Hl- 

^ . The Mathematical Model for the Missile Intercept Problem 

In this section the mathematical model for the missile 
Intercept problem solved in Section V is derived. An air- 
to-air missile launched from an attacking aircraft must 
intercept a constant-velocity target. The missile is 
restricted to maneuver in a plane. The orientation of this 
maneuver plane in three dimensional space is defined at 
launch as the plane containing the position of the missile 
at launch, the position of the target at launch, and the 
velocity vector of the target. 

The assumptions applied to the problem include those 
presented in Section A.l plus the following: 

a. the initial velocity vector of the missile lies in 
the maneuver plane, 

b. the attacking aircraft is tracking the target at 
missile launch so that the missile's initial velocity points 
at the target at t = 0, 

c. the target moves with constant velocity, 

d. components of out of plane forces perpendicular to 
the maneuver plane are ignored. 
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Figure 40 presents a view of the problem in the maneuver 



plane . 



TGT 




Figure 

Missile Coordinate System 
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Four axis systems are drawn in Figure 40. They are; 



a. (X',y) 

b. (X ,Y) 

c. (x",y") 

d. (x',y') 



a fixed inertial axis system in the maneuver 
plane with the origin at the missile launch 
point; 

a Newtonian reference system in the maneuver 
plane v/ith the origin at the missile launch 
point at t = 0; after launch the origin 
remains fixed with respect to the target 
(it moves with velocity v„ with respect to 
the X'Y' system); 

a non-rotating axis system fixed to the 
missile center of mass; 

a rotating axis system fixed to the missile 
center of mass; the x' axis is oriented in 
the direction of the missile's velocity 
vector. 



The systems X'Y', XY, and x"y" are oriented in the maneuver 
plane so that the axes O'X', OX, and cx" form the intersection 
of the maneuver plane and a horizontal plane. These axes are 
chosen so that the target's initial X, X', and x" positions 
are positive. The axes O'Y', OY, and cy" are chosen so that 
the conponent of missile weight inthe maneuver plane is acting 
in the negative Y', Y, or y" direction. All angles are 
positive as they are shown in Figure 40 in the counter- 
clockwise direction. The remaining notation is: 



a. forces; 

N = normal aerodynamic force , 

A = axial aerodynamic force, 

T = thrust, 

Th = normalized thrust, 

W = component of missile weight in the maneuver 
^ plane, 

b. angles; 

a = angle of attack ^ 

6 = missile flight path angle, 

Y = target flight path angle. 
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c. other quantities; 



V = missile velocity, 

= target velocity, 

M = missile Mach number, 

= target Mach number, 

Wg = missile gross weight, 

Cjj = normal force coefficient, 

= axial force coefficient, 

(j) = angular velocity of the x’y’ system with 

respect to the x*'y" system, 

m = missile mass 

> 

g = gravitational constant, 

S = missile wing area^ 
p = air density, 
n = load factor, 
a = speed of sound. 

The equations of motion are derived following the methods 
used in Section A.l. Equations (A.l) thru (A. 5) are 
identical. The angular velocity of the x’y' system v;ith 
respect to the non-rotating frame x"y" is given by 




The product 



0 ) X V = V e e„, 



(A. 65) 



(A. 66) 



is formed. Summing forces in the x' and y' directions, we 
obtain from equation (A.l) 



I \ 
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« 




(A. 67) 



Tcosa - Nslna - Acosa - W sin© = mv , 

c ’ 

Tslna + Ncosa - Aslna - W^cos© = mv © . (A. 68) 

Axial and normal force coefficients are defined by the 
expressions 



A = Ca I pv^S , (A. 69) 

and 

N = Cj^ I pv^S . (A. 70) 

Substituting equations (A. 21), (A. 22), and (A. 23) into 
equations (A. 69) and (A. 70) and transforming the results 
Into state variable format, the state equations become 



M = 



aW 



Thcosa - 



KpSa 

2W 

e 



2 

M C^cosa 



- - 

e 



gw 

aW. 



sin© 



(A. 71) 



• ^ gjjsjjiq _ ggSa . geSa 

® aW M 2W « 2W 

e e e 



MCj^cosa - 



gW 



c cos© 
aW M 



(A. 72) 



Two additional state equations are required to impose end 
conditions on the relative positions of the missile and 
target in the optimization procedure. They are 



X = aMcos© - aM^cosy 



(A. 73) 



Y = aMsin© - aM^siny 



(A. 7^) 



•' \ 
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The states are missile Mach number M and missile flight 
path angle 0. The control is missile angle of attack a. 
Normalized thrust is given by 



Th = 1 


v| 

-p 




(A. 75) 


Th = 0 


t > 






where t_ represents engine 

D 


burnout . 


The 


follov;ing inequality 


constraints are Imposed; 




• 




v| 


® < 0... 
— M 




(A. 77) 


"m - g 


(6M) < 




(A. 78) 



Equation (A. 78) represents a structural load factor limit. 

The axial and normal force coefficients are given as 
tabular functions of Mach number and angle of attack. 
Parameters considered constant for the problem are the 
gravitational constant g, maximum thrust Thj^, the speed of 
sound a, missile weight W^, air density p, missile wing area 
S, the Mach number of the target M^, and the flight path 
angle of the target y. 

In order to properly define the problem, it is necessary 
to perform several manipulations in analytic geometry. 

First, the three dimensional positions of the missile and 
target must be specified at launch. Second, the velocity 



< \ 
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vector of the target and the Mach number of the missile at 
launch must be specified. Once this is done it is necessary 
to : 



•a.' identify the maneuver plane, 

b. identify the XY coordinate system, 

c. calculate the target coordinates in that system, 

d. calculate the target flight path angle Yj the initial 

missile flight path angle 0(0), and the component of the 

missile weight acting in the maneuver plane V/ . The 

c 

optimization procedure can then be commenced. 

To accomplish these calculations an initial coordinate 
system is established in which the problem can be easily 
visualized. The origin is situated at the missile. The OX 
axis is positioned in the horizontal plane. The OY axis is 
positioned in the horizontal plane such that the target has 
no Y‘ coordinate . The OZ axis is positioned in the vertical 
plane such that a target which has an altitude advantage 
over the missile has a positive Z component. This coordinate 
system is shown in figures ^1 and ^2. The angles 6^ and 3^ 
are defined as shown above. The following relations may be 
written: 



a = R^i + h^k 



(A. 79) 




(A. 80) 



tan6„ = — 
T m^, 



(A. 81) 



l6l 



z 




measured in the 
vertical plane 
containing 



Figure 



M TGT 




Figure ^2 

Initial Missile Coordinate System 
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m 



tan3^ 



y 



(A. 82) 



m 



X 



With the problem defined in the coordinate system 
shown in Figures and H2 it is now necessary to transfer 
the problem to the coordinate system used in the optimization 
procedure. That is, it is necessary to identify the maneuver 
plane and the XY coordinate system. To this end, a vector 
normal to the maneuver plane is 



To establish the X axis a vector is required v;hich is in 
both the maneuver plane and a horizontal plane. Such a 
vector is 



To establish the Y axis a vector is required which is in 
the maneuver plane and perpendicular to the X axis. Such 



N = a X M, 



(A. 83) 




X = N X k 



(A. 85) 



(hTHix.-RTinzt )1 + hipm^.j . 



(A. 86) 



a vector is 




\ 



i 



y = N X X 



“ —Riph,piTij^ , 1 + R^niy , ( h,j|in^ , R,pni^ , ) J 

=[h^^my,^+(h^m^,-R^m^, )^]k . 



(A. 87) 
(A. 88) 



The angle {() between the maneuver plane and a horizontal 
plane is required and is given by 



COS<|) 



N 


• k 


|N1 


|k| 







(A. 89) 






[(h,pm , )^ + (h^m^,-R^m^, )^ + (R^m , )^]^ 



(A. 90) 



The missile weight component in the maneuver plane W may 

c 

be found from 



Wq = WgSinct. • 



(A. 91) 



This is Shown graphically in Figure ^3. 




maneuver plane 



Figure ^ 3 

Missile Weight Component 
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The initial target coordinates are 



X^(0) = 



PROJ^a I 









(A. 93) 



Y^(0) 



iPROJ^a 



(A. 95) 



= + 









The sign of Y^(0) is resolved by: 



if > 0 , 


Y^ is 


positive ; 


if h^ < 0 , 


Y^ is 


negative . 



The initial missile flight path angle 0(0) is 



tane (0) 



Y^CO) 

34X07 



(A. 96 ) 



by assumption b at the beginning of this section. Before 
proceeding it is necessary to insure that the X and Y 
vectors given in equations (A. 86) and (A. 88) have the 



correct sense. This may be checked by observing the sign 
of PROJj^a which should be positive and the sign of PROJya 
which should be positive if h,p > 0 or negative if h,p < 0. 
After the senses of these vectors have been checked and 
altered as required, the target flight path angle y 
calculated by 



5r*x 

cosy = — . 

~ I 

The possible range of y is 



(A. 97) 



-TT £ y £ TT . (A. 98 ) 

If cosy is positive, then 

- I < y < I . (A. 99) 

If cosy is negative, then 

^£y£TT or • (A. 100) 

To find which inequality applies in (A. 100) the quantity 



= ~T*^ 

^ ' I??tII!I 

is formed. If k is positive, then 

0 < y < n 



(A. 101) 



(A. 102) 
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If k is negative, then 



TT < Y < 0 



(A. 103) 



This logic completes the set up of the problem In the 
maneuver plane. 



APPENDIX B 



TABULAR FUNCTIONS 



In this Appendix the tabular functions used In the 
problems are presented. 

1 . Three Dimensional Plots 

Three dimensional plots of all tabular functions are 
presented here. Figure gives the lift coefficient as 
a function of Mach number M and angle of attack a for the 
the supersonic fighter aircraft used In the aircraft 
problems. Figure 45 gives the drag coefficient as a 
function of Mach number M and angle of attack a for the same 
fighter. Figure 46 gives normalized maximum thrust Th as a 
function of Mach number M and altitude H for the supersonic 
aircraft performing the minimum-time climb In the problem In 
Section VI. Figure 47 gives the normal force coefficient Cj^ 
as a function of Mach number M and angle of attack a for the 
air-to-air missile used In the problem In Section V. Figure 
48 gives the axial force coefficient for this missile. 

2. Tables 

Following each plot a condensed version of the table used 
in the computation Is presented. 



\ 
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Figure 
= f(M,a) 
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TABLE 11 

LIFT COEFFICIENT FOR A SUPERSONIC AIRCRAFT 



O OOOIAOIA NTOrOOiAO— ^OOGOOmOGOO— ‘O— 'OrOO rnoOOsOOOO — lO 

O vOodoosJ4-o orsj^nmr-ocsjnr-icooa'vOor^rsjirMno-Hoo r-iococDOu^or^r^O 

CO oo^coin^vo vO<Mco>ror-mo>iAOvOomocor«-o^vOcoor-iooo>ro>oO'rsy GOO 
• (^JrHr-^^^ooo-HrH<M^o^o>^’^^• vOf^r^ooGOaDo>o>0'0>oo OoO'O'ooGOh-vO LfMn 

o 

OOOOOOOOOOOOOOOOOOOOOOOOOOO— irH-HOoOOOOoOO 

I I I I t I 



O OcoOrOO'T ChO^O 'ror^0rn0—i0rnor00--<00>0>ro vOOC>OoOChO 'I'O 

o nO'OQoocnjO' 0 '*-‘^<Mr 0 inv 0 r^> 4 - 0 ^pHO>rv 0 i^' 0 — •(M^ncooo r^ocoorncoNr 

r- oo^r^iArHin Lrsrsjco'roNOf\joo>rcMnO'Tco(\»sOor<uAr^orO'-^o^u^fNjr^— ‘incoo— < 

• cNj^^'-<^ooo«-ir-i(\ifr>rn >TN^Lrunvor^r^r^Goooo'0'o>^oo oa^cr^g^cocor^vO'^in 

o • 

OOOOOOoOOOOOOOOOOOOOOOoOOOOr-HrH rHoOOoOOO oo 

I I I I I I 



o OrOOvOOm >t*OOOoONOO^OrOOmo<J'OOOrHOvOOO>OC\lOC 7 'OcooOo 

o NOvOr^coOuo c^ooooooO'OOf'jfNjcOfMCNjomno— iirvcoo<^oQ'c\j'4*ocoofn GO 

vO mrsjoc>r coco csjvOO^c^omr^r-im mo coco rOvOOOsim 

• osjr-i'— <•— ‘•"HOoOr-HrHOsjrnro 'Tv^-iAiAvO'Or^r^Gocoooa'O'O^oooocT'O'CDQor^r^vOin 

o • • 

0 oo oo O oo O oo O O O O o O O O oo O O O O O O r-i r-i pHrHO o o oo o o o 

1 t M I I 







u 


o 


< 


o 


H- 


lA 


CCh^ 


• 


IXJ< 


o 


CO 




XLL 




Oo 








uu 


o 


X_j 


o 


oo 


M- 




• 


2 :< 


o 


n 11 




aicc 




UJUJ 


o 


h-K 


o 


LULJU 


CA 


2 TX 


• 


<< 


o 


cCcc 




<< 




QlO. 




o 


o 


< 


o 


I--J 


CM 




• 


oo 


o 






HHJ— 




c^otl 




OLXJ 




x> 


o 



o 

.-I 

o 



0^^0^00^0 OOOOOOO'-HOrOOtnO'TOOOr-iO OO0>O<MOQ'O00OOO 
vO^f^COOCO 0'OOOOOOOvOOoj(\loOf\J(MOLncOr-iin:OOrno^ CNj>r O CO O rn ro 
incsjcDvrOvOrNiconia' Nr\-r'rooy<N»'>ou'rsjLn:^.-HLr\ mo'^>f’'iomooosjm 

f vj ^ rH ^ ^ o o o ^ rH ( sj rn r<^ ^ lA m sO >0 h* a:^ CO CO a> o”' o o o o tr 30 CO r- vO in 



OOOoOOOOOOOOOOOOOOOOOOOOOOOrH—j^^OOOOoOOO 
\ \ \ \ \ \ 



omOvoom ' 4 *ooooo>oo-Horoocno>rooor-io ooooosiooocoooo 
vOvOf^cooco O'OOooooO'OoosirsjTOOsirsJomoo— lu‘^coo^’Oo^(^J^rocoo mco 
oCM^>t-^iA lAosjcOvroocNJooro J'moo(^J^oc7^(^J’A^-^-lln^">o^^ocomoocM^o 
<M^ rH ^ r-i o o o ^ --H ro 'I' lA LA vO o CO 00 CO CA O' o o o o O O' 00 CO r- h* o in 



OOOOOOOOoOOOOOOOOOOOOOOOOOOr-irHrHrHOOOOOOOO 

I I I I I I 



omoNoOrn >i'ooooo>oo-^orooiAO>^oa>o-HO oo O'OcoochOQoooo 

vO 'Oh- GOO CO OOOOOOOOOO AJIMCOAJAJOLACO— ‘lAcOO ^0OC^ C>sJ >4" O T> O CO 
OOf^<t**-»LA 'A ^MOO>^0'Of^Jco^OC^'4'^J^^OOO(^J>0'7^(\|u"^f^^LAfOOvO^OCOo> oOosjm 
CNJ r-l r-< ^ O O O '-ICNJ rO ro st LA LA >0 nO h- 00 CO 30 v/' O O O O C> c* CO 00 h* h- nD LA 



OOOOOoOOOOOOOOOOOOOOOOOOOOOrH^r-irHOOOOOOOO 

I I I I t I 



OrOOvOOrA >tOOOOOvGO^Or0OlAO4‘OOO'-<O0OChOC\JOChOG0OOO 
vO'ON'OOOOO ChOOOOOOO oo AJOslOOCvJAgOLAOO— tlAOOO rOOOCsJ'J-OoOOrAGO 
O lA<\JcO>^0'0(NJOOrOC>'1’OrACOOsjOOc\|lA'^^lArOOvOrAOOrAvOO(NjrA 

CM r -4 rH rH —I o o o —I •-« iNj ro ro >t 'I' LA lA vO vO CO 00 CO c^ o o o o o (j> ch 00 CO vO LA 



OOOOOOOOOOOOOOOOOOOOOOOOOOO^— «rHrHOOOOOOOO 

t 1 I I t » 



OrOO OOrO M-OOoOO>0O^OfAOiAOstOC^Or-HO>0O GhOCMOO'CCDO oO 
vOOr^ooooo a^ooooooo^ooc^JAJcofMf^JOu^cor-^l^^coo ^ooc^(^Jvtocoo^A *o 
OOr- M'^iA LAC^JOO^^-OvO(^JCOfAC^'4'C^^ACO.■^J^UC^c^JLA^••-HJ^^OOvOxA 33 oA >0 O CM rA 
C\jr-<r-^.-<r-^oOO'-<^rMrArAM" >J*lALAvO or^^-COOO COO'* OO^OO OO O^U'COcO^-r^ 



OOOOOOOOOOOOOO^OOOOOOOOOOOO— «— «^rHOOOOOOOO 

I I I I I I . 



OrAO OOrA vtOOOOOvOO^OrAOLAON^OO'O— <OvOOO>OcMOO'OQOOOO 
vOvOr^COOCO O'OOOOOOONOOCMCMOOCMvMOLAOOr-iLnOOO rAOO>CMM-OCOO ^'JO 

o . ocrr^M'^iA lACMOOM-OvUCNjcorAu^'ro' A>cDvM^uocMuv^r-iLno^ONartiajrANOOCMrA 

• CMrHrH^— •OOOr-t-HrsjrArAM'M’LAlA QO 33 COU' ^J>C^OO OOO^ O' OJcOl^r- nQLA 

o 

O OOOOOOOOOOOOOOOOOOOoOOOOOO*-HrH ^^OOOOOOOO 

I I I I I I 



OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 



nO LA rA (M — I o — * CM rA Ln vO A 3 O' o — ICM rA 4 * ir^ nO 00 0 > O — I CM rA M" LA vO A 3 O' O — I ^ 

I I t I I I ^r-ir-|^^^-H^rHr-iCM(M(MfMCMCMCMCM<MCMrACArA 



\ 



\ 



170 



TABLE 11 (continued) 

LIFT COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 11 (continued) 

LIFT COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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Figure 

Cp = f(M,a) 
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TABLE 12 

DRAG COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 12 (Continued) 

DRAG COEFFICIENT FOR A SUPERSONIC AIRCRAFT 



O s0 0^rT|cr>'rc\jc\ic\J4*00frjOs0inmr^0'^0'v0N}'^Lnr^^'Oc\jO0^Or-^incMnc\i^r-^c\jm 

h- sOiAin't>J" 4 - 4 ''tNt^iALnor^ooo'Oc\iroLnh- a^^rnvooo^>j-NOorooa^f<^r^^Lna^cri 

• OOOOOOOOOOOOOOOO^^^^^-HC\|C\Jc\jC\Jrnr 0 r 0 nO^- 



OOOOOOOOOOOOOOOOOOOOOO^OOOOOOOOOOOOOOOO 



o o vOh-ONOrnc\imvoor^ ^o^oc^f^c^^^^^c^foc^ >o 'Oco^ vOmro4' sO'-^ooh-r^o ^^’Oc^C7' 
vO voinm'TN}'>T'TNi-LninvOh-ooo^rninr^oc\JLn^^mcof\jvOQ moo roo>ro 
t oooooooooooooo^'-'»^*^^f\ifNJi'Of\immm<f >n^\L'v^\^£)^o^-^-ma‘C^o 



oooooooooooooooooooooooooooooooooooooo^ 



o c^^-f-^r-^^-ooLnoo^-•-^r-^^-o^^oc^commc^^c^c\l•-^^o^^m^m^-^-f\lmC'f\lo^^“^ 0'0 

o oo* 4 "m'tr^mc\imh-' 4 "m^oom>^-o^ 2 Dr^o* 4 'rsic\jN^'C'f^r^OLr>m> 3 ’r^rj^cvLn^c^^ 

m ooh-vOLH^t ' 4 ‘^^>J"i^o^-oooc\J 4 ’r^^c\i'Oa'mf^^moin^vOc\ioo>r^aotnrvjor-o 

• ooooooooooooo^'-*^^^fNjfNJf^Jmm>r >^’ir^msOvor^r^oocr'a'o^fNjr\jm 



0000000000000000000000000000000000 ^^^^-^ 







o 


O 


< 


o 


h- 


4- 


(XH- 


« 


m<i 




CD 




SILL 




OO 




2 




m 


o 


xm 


o 


oo 


m 




« 


5I< 


rH 


H 11 




CCOC 




miu 


O 


h-H- 


o 


mm 


CM 


XX 


• 


<i< 


rH 


<X(X 




<< 




CL CL 





C 7 'vOin:jOiruAooLnincoinvoa' ^r^r-iooa>mo^ vom^cM^oomr-imooooom^mcM^o 
CM vt O' h- (jN c\j nJ- o r- o> 'I- CM 4 * O' 00 O' m lA o 00 o tn m m o ‘ 7 ^ ^ vO > 3 " vO o f\j Lf^ c\J 

ooovOiA>i->r^rM-<r»AvOcoo^M'h-o>roor‘jr^r- 4 r-(Mcj -r^i^iACMOoor^'OinM“>r 
r-i o o o o o o o o o o o r-i r-» ^ ^ CM CM CM m m 'T 'T min no r- 00 O' o o ^ 'M m M" cA vO r- 



OOOOoOOOOOOOOOOOOOOOOOOOOOOOOO^i-<^f-<»HrHf-Hr-|r-i 



CM CM O' h- m m m o^ CM CM 00 1 A m CM 'H CM > 3 - vO o o o m ^ o o m o o o > 3 - <M CM 

I — 3 ‘v 0 '-*^mf 0 m ^^'0 4 'h“m' 3 -a'oo^co 0 ' > 3 -' 3 "^-m> 0 (\jcM-OM"'OCMcM^-mco> 3 "mO T' 
r-i h- vO m ^ ^3- ^ m no O' r-i h- o > 3 ’ O' m 00 nJ- o 'O m o CO vO nJ- m CM CM CM CM m > 3 " vO CO m 
r-* o o o o o o o o o o o r-i ^ r-i CM c Ni CM m m > 3 - m m nO 00 O' o CM m > 3 ' m vO 00 o ^ 



OOOOOOOOOOOOOOOOOOOOOOOOOOOO'-C^'-*^— ' ^r-l' vJCM 



v 0 h-^ 0 ' 00 v 000 o>> 3 -h- vOm OOO' 00 ^h-co 0 '-»rOCMO'm^moOr-iCMO'M’vOiA^mvO' 3 - 
O' m CM m s3“ Nt h- 4 * m CM m Qv o m m o o Nt* -o O' r*- O' r- o "0 h- M' m o ^ vO o ^ o ' 3 - m 

CM O 00 vO m 4 " 4 ' lA nO CO O CM ' O O m CO m 00 4 * O 4 - CM O O' CO OO CO O' O CM ' 1 * o nO 

r-Cf-iooooooooo*-^ .-<r-i^cMcMmm 4 'm'-n>or^crocoo' 0 »^cM 4 'mNOh-coo— *on 4 - 



000000000000000000000000000 ^^'-*'-*r-i^r-i^CMCMCMCM 



-» o 

<t o 



z< • 
Oo 

isjl— t 

o^q: 

am 

x> o 
o 
o 



O' o lA r-i o ^ lA r-i o r-* m o o> o m O ' 00 i-i vO 4- <r h- m o o m 00 m m 00 CM o O' ^ vO m CM 4" 

COO'D CO mh- 4‘h- AoovOOoocAcM'OOrHCMi^oo'3'mcM4“r-imom^4'mh-moo'4-4’0' 

m r-ico vo m 4* 4 - 4- m vo 00 ^ m ^m o o CM CO Ai m O' O' O' o r-i m mco ^ m o 4" o 'D CM 
•-H r-i o o o o o o o o o r-i rH CM CM m m 4" '3* m vO c-o 00 O' o CM m 4" m o 00 O' ^ CM 4" A r- 



O OOO OOO O O O O O OOo O OOOO oo OOO Or-Ir-I ,-1 CMCM CMCMCM 



4 " CO A h- CM r-i 4 - ^ CM m CO 4 “ m O' O'* m m O' vO 00 m m 'O m 4 - O ' 00 •-< CO CO m ^ m O' o 4 " 
m m r-icMO' '-coo ^ O' CM m moo 'O 4 " o m >0 4 " cor^ 'D nocm 4 * 0 cMO m ^ 4 * 4 * 00 CO 

m o 00 o 4* 4- m 4* 4* sO CO o m ^ A o sO CM O' 'O 4* CM rH »-H ^ r-< CM 4 - nO O' CM o o m o 'O Aj 

rHr- 4 OoOoooooO'- 4 »-^«-HcMCMf 0 m 3 - 4 -mvOh-' 00 'o*-<cMm 4 -mr^coo^f'") 4 ' or^ 



OO oo o oo o oo Oo oo OOO OOo o oo oo 'H CM CM CM CM CMCM 



o cMO'C^l-l'OCMr^cM^O'-lO'a'CM^- 4 ■mmc^mm 4*h-cMOO'CM'OCMrHmsOCMoomr^ 4-4'm 

o r- m 4 ‘ 00 ' rO'-iro QN o 4" m 4 - nocm CM Dm CO m^cMCMm >3- Num '3' O' ^cMo CM u) com cMm 

O' a'h-u\ 4 'CMCMcMfMcM>rmh- LT'CMmu'mr'-cMr^mo' DmocovL>m>Trr>m 4 ’mDr- J'cmuaoo 

• ooooooocDooooo^r-i^cMCMmm>T 4 'mDr-r-Loa' 0 ^cMm>r mD^-CT'O^ 

O 

0000000000000000000000000000-^^^--»'-*--*^^--*CMCM 



OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 

sO m 4" m CM r-i o CM m 4" m D 00 O' o CM m 4* m D 00 O' o r-* CM m 4" m vO 00 O' o r-i CM 

I I I I I ). r-l^r-lrHr-lr-|rH^^^CMCMCMCMCMCMCMCM‘MCMmmm 



\ 



\ 



175 



TABLE 12 (Continued) 

DRAG COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 13 (continued) 

MAXIMUM THRUST FOR A SUPERSONIC AIRCRAFT 
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TABLE 13 (continued) 

MAXIMUM THRUST FOR A SUPERSONIC AIRCRAFT 
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TABLE 14 

NORMAL COEFFICIENT FOR AN AIR-TO-AIR MISSILE 
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TABLE 14 (continued) 
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TABLE 15 (continued) 
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APPENDIX C 



INTERPOLATION 

In the optimization problems solved herein the aero- 
dynamic data is given in tabular form. The dependent 
variable D is given as a function of two independent 
variables M and a in all cases. Excerpts from the tables 
used in computation are presented in Appendix B. 

2 2 

For a given M and a quantities D, ^ ^ j , 

3 2d 

and are required by the optimization algorithm. 

Parabolic interpolation is used to obtain these quantities. 
In this Appendix parabolic interpolation for two independent 
variables is derived. 

1. Parabolic Interpolation in the Plane 

To apply parabolic interpolation to a tabular function 
of two Independent variables the nearest point in the 
tables to the given point (M,a) must first be found. It is 
assumed that the tabular data is given at constant intervals 
AM and Aa in the independent variables. The nearest point 
given in the tables and the surrounding eight points are 
required in the interpolation and are shown in Figure 49. 

The parameters 0 and 4> locate the point (M,a) from the 
nearest tabular point (M ,a ). If (M ,a ) is the nearest 

3 3 3 3 

point then 

- J 1 1 J (C.l) 
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Figure 49 



and 

" I 1 0 1 I (C.2) 

The inequalities (C.l) and (C.2) hold unless the nearest 

point (M ,a ) is on a border of the table. In this case 
s s 

(M ,a ) is chosen one point in from the border. The 
s s 

parameters 0 and then satisfy 

- 1 < 4) < 1 (C.3) 



and 



1 < 0 < 1 



(C.4) 



V/rltlng Taylor series expansions Including up to second 



order terms for the eight points surrounding (M ,a ) , v;e 

s s 

obtain 



D(Mg+AM,ag+Aa) - ) = D(Mg,ag) 

3D 



(C.5) 



+ am 
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ID 


+ 


(AM)^ 


3^D 


+ 


(aci)2 


3^D 


3a 


s 


2 


3M^ 


s 


2 


3a^ 




“s+1^ 


~ D(M ,a 
s 


s> 
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3M^ 


1 
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+ AMAa 



3^D 

3M3a 



(C.6) 



+ AMAa 



3^D 

3 m 3 a 



(C.7) 



- AMAa 



3^p 

3'M3a 



(C.8) 



- AMAa 



3M3a 



D(M^+AM,a^) = D(M ,a ) = D(M^,a^) + AM 

o o S'J.S 00 



s " s 



3D 



. (AM)^ 3^D 
2 . 3M^ 



(C.9) 
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s s 






D(M ,a_ +Aa) 

o o 






^ I? 



V(M^,<x^-La) = DCM^.a^.^) = DCM^.a^) - ia ff 



(AM)^ 


3^D 




2 


3M^ 


s 




(C.IO) 


(Aa)2 


2 

3^D 




2 


2 






9a 


s 




(C.ll) 


(Aa)^ 


9^D 




2 


9a2 


s 



(C.12) 



Subtracting equation (C.IO) from equation (C.9), we obtain 



3D 



- DCM^.i.a^) = 2AH ^ 



or 



9D 

3M 



2AM 



(C.13) 



Subtracting equation (C.12) from equation (C.ll), we obtain 



= 2''“ II 



or 



ID 

da 






2Aa 



(C.li4) 
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Adding equation (C.IO) to equation (C.9), we obtain 



^ = 2D(H^,a,) + UM)^ ^ 



or 



a^D 



9M 



(AM)^ 



(C.15) 



Adding equation (C.12) to equation (C.ll), we obtain 



D'«s>“s-H> ^ D(M^.a3.i) = SDCM^.a^) + Ua)^ ^ 

9a 



or 



a^D 

aa^ 



(Aa)^ 



(C.16) 



Subtracting equation (C.7) plus equation (C.8) from equation 
(C.5) plus equation (C.6), we obtain 






s-l 



“s-1^ 









^AMAa 



2 

a D 



aMaa 



(C.17) 



or 



2 

a D 



aMaa 



lliMAo 




i 



A Taylor series Is v/rltten for the point (M,a). Using 
equations (C.13), C.l4), (C.15), (C.l6), and (C.17), we may 
reduce this series to 



D(Mg+6AM,ag+<()Aa) = D(M,a) = D(Mg,ag) + | - D(M^_^,ag)] 









(C.18) 






Rearranging, we have 



D(M,o) = ^ 



(C.19) 



+ (l-e^-<J)^) D(Mg,a^) 
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Equation (C.19) is the expression used to interpolate for 
the value of D in terms of the surrounding none tabular 
points. To obtain expressions for the required partial 
derivatives, observe that 



M = M + 6AM 
s 



(C.20) 



and 



a = a + (bAa 
s 



(C.21) 



Therefore , 



dM AM 



(C.22) 



and 



= Ji. 
da Aa 



(C.23) 



The chain rule for partial derivatives yields 



(Mg+6AM,ag+<bAa) 






aD(M,a) ae 

9 6 9M 



(C.24) 



Taking the partial derivative of equation (C.19) with 
respect to 6, we obtain 
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- 2e D(Mg,ag)] . 

Using similar procedures, v;e may derive the remaining 
expressions , 



3D 

3a 



(M.a) > ^ [ f ^ 

- 2* D(Mg,Og)] 



■ f 

(C-26’ 



^ (M,o) = 



1 

(AM)^ 






(C.27) 



2 

^ (M,a) = CD(Mg,ag_^) - 2D(Mg,ag) + D(Mg,ag^.;^)] (C.28) 

3a (Aa) 



3^D 

3M3a 



= W C°(«s-rVi> - ^ 



(C.29) 
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APPENDIX D 



EMPIRICAL RELATIONS 

In the problem treated in Section VI empirical relations 
are used for 



a = f(H), 



(D.l) 



a = f(H), 



and 



M^ = f(H) 



(D.2) 

(D.3) 



where the parameters are air density ratio (o), speed of 
sound (a), maximum Mach number and altitude (H) . The 

air density ratio is defined as 



A 



o = 



Po 






where a is sea' level standard day density. This Appendix 
presents these empirical relations and compares the values 
obtained from these relations xvith standard atmospheric 
conditions . 

1 . Air Density Ratio 

The empirical relation used for air density ratio is 

-c-|h -c«h 

a = e +c^h (D.5) 



; 
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where 




and 

= 36,089 ft. 
= 1,5^100 
C 2 = 1.804ii5 
= o.iaso 



(D.6) 



(D.7) 

(D.8) 

(D.9) 

(D.IO) 



Figure 50 is a plot of the values of a obtained from 
equation (D.5) compared to those obtained from standard 
atmospheric tables. 

2 . Speed of Sound 

The empirical relation used for the speed of sound is 



a = a^(l-c.^h) , h < 1 (D.ll) 

= 971 ft. /sec. , h > 1 (D.12) 

where the parameter a^ is the speed of sound at sea level 
on a standard day; that is 

a^ = 1116.89 ft. /sec. (D.13) 

and 

= 0.1331 (D.14) 

Figure 51 is a plot of a vs. H. The expressions (D.ll) and 
(D.12) duplicate exactly the values obtained from standard 
atmospheric tables. 
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air density (slugs/ft. 



2.4 




Figure 50 

Air density vs. altitude 
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Figure 51 

Speed of sound vs. altitude 
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3. Maximum Mach Number 



The empirical relation used for the maximum Mach 
number of the aircraft for the problem solved in Section VI 
is 



Mj^ = 2.1 - l.le 



-2.Hh 



(D.15) 



Figure 52 is a plot of equation (E.9) along with the actual 
restrictions of the aircraft under consideration. 




Figure 52 

Placard Mach number vs. altitude 
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APPENDIX E 



A CONVEXITY THEOREM 



In this Appendix the follov/ing theorem on convexity Is 
proved. 

Theorem 1; If f(x) Is convex on v/here x e R^, and f ^ 0, 

K n 

then f (x) Is convex on R v/here K Is any positive Integer. 

This theorem Is proved by mathematical Induction. First, 
the follovflng theorem Is proved. 

Theorem 2; If f(x) Is convex on R*^ v/here x e R^, and f ^ 0, 

2 n 

then f^(x) Is convex on R^. 

Proof : The function f(x) Is convex If 

f[XX 2 + (l-X)x^] £ Xf(x 2 ) + (l-X)f(x^) (E.l) 

for all x^, X 2 and X e [0,1]. Squaring both sides of 
Inequality (E.l) we have 

f^[XX 2 + (l-X)x^] £ [Xf(x 2 )+ (l-X)f(x^)]^ . (E.2) 

The sense of the Inequality (E.2) Is retained as f ^ 0 . To 

2 

prove that f (x) Is convex, It must be shovm that 

f^[Xx 2 + (l-X)x^] < Xf^(x 2 ) + (l-X)f^(x^) . (E.3) 

Observing Inequality (E.2), (E.3) Is seen to be a true 
Inequality If It can be shown that 
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[Xf(x 2 ) + (l-X)f(x^)]^ < Xf^(x 2 ) + (l-X)f^(x^) . (E.iO 

To show that Inequality (E.4) is true, we proceed as f ollov/s . 
Since X e [0,1], it is true that 

X[f(x2) - _< [f(x2)- f(x^)]^ . (E.5) 

Expanding inequality (E.5), we have 

X[f(x2 - f(x^)]^ 1 .(E.6) 

Rearranging inequality (E.6), we have 

X[f(x2) - f(x^)]^ + 2f(x2)f(x^) - f^(x^) £ f^(x2).(E.7) 
2 

Subtracting f (x^) from both sides, we obtain 

X[f(x2) - f(x^)]^ + 2f(x^)[f(x2) - f‘(X]^)] _< ^*^(^2) - • 

(E.8) 

Multiplying inequality (E.8) by X , we have 

X^f^(x2) - 2X^f (X2)f(x^) + X^f^(x^) + 2Xf(x^)f(x2) (E.9) 

- 2Xf^(x^) _< Xf^(x2) - >^f^(x^) . 
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Adding to both sides of inequality (E.9), we obtain 



X^f^(x2) - 2X^f(x2)f(x^) + 2Xf(xj^)f(x2) + 

(E.IO) 

-2Xf^(x3^) + X^f^(x^) < Xf^(x2> - + f^(x^) . 

Simplifying inequality (E.IO), vie obtain 

CXf(x 2 ) + (l-X)f(x^)]^ < Xf^(x 2 )+ (l-X)f^(x^) . (E.ll) 

This is the inequality we set out to shov/. The theorem is 
proved. 

Second, the following theorem is proved. 

Theorem 3* If f^(x) is convex on where x e R*^ , and f > 0, 
then f^^^(x) is convex on R^. 

Proof; It has already been shown that if f(x) is convex and 
f ^ 0, then f (x) is convex. It is now assumed that 

f^[Xx 2 + (l-X)x^] < Xf^(x 2 ) + (l-X)f^(X 3 ^) . (E.12) 

Multiplying both sides of inequality (E.12) by f[Xx 2 + (l-X)x^] 
a positive quantity, we obtain 

f^'^^[Xx 2 + (l-X)x^] < [Xf^(x 2 ) + (l-X)f^(x^)]{f[Xx 2 +(l-X)x^]} . 

(E.13) 
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Substituting the expression (E.l) into inequality (E.13), 
we have 

f^‘^^[Xx2+(l-X)x^] < [Xf^(x2) + (l-X)f^(x^)][Xf(x2) + (l-X)f(x^)]. 

(E.liO 

K+1 

To prove that f (x) is convex it must be shovm that 

f^‘‘'^[Xx2+(l-X)Xi] < Xf^'^^(x2) + (l-X)f^'^^(x^) . (E.15) 

Observing inequality (E.l4), (E.15) is seen to be a true 
inequality if it can be shown that 



[Xf^(x2)+(l-X)f^(x^)][Xf(x2)+(l-X)f(x^)] < Xf^'^^(x2) + (l-X)f^'^^(x^) . 

(E.16) 

To show that inequality (E.l6) is true we proceed as follows. 
The expression 

Cf^(x2) - f^(x^)][f(x2) - . (E.17) 

is always a positive number because the signs of the expres- 
sions in parentheses in expression (E.17) must be the same. 
The following inequality holds 

X[f^(x2)-f^(x^)][f(x2)-f(x^)] < Cf^(x2)-f^(x^)][f(x2)-f(x^)]. 

(E.18) 



> \ 
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Ejcpanding and multiplying by X, we obtain 






X[f^'^^(x2)-f^(x^)f(x2)-f(x^)f^(x2)+f^'^^(x^)] 



(E.19) 



K4* 1 

Rearranging inequality (E.19) and subtracting 2Xf (x^^) 
from both sides, we obtain 



X^[f^'^^(x2)-f^(x^)f(x2)-f(x^)f^(x2)+f^‘*'^(x^)] 



(E.20) 



+ X[f^(x^)f(x2)+f(x^)f^(x2)-2f^'*'^(x^)] < X[f^'^^(x2)-f^'^^(x^)] 



Further expansion and rearranging yields 



X^f^'^^(x2)+Xf^(x^)f(x2)-X^f^(xj^)f(x2)+Xf(x^)f^(x2)-X^f(x^)f^(x2) 



- 2Xf^"''^(x, )+X^f^‘‘‘^(x, ) < Xf^"''^(xp)-Xf^'''^(x, ) . 



(E.21) 



K+1 

Adding f (x^) to both sides and rearranging further, we 
have 



X^f^'^^(x2)+X(l-X)f^(x^)f(x2)+X(l-X)f(x^)f^(x2) + (1-X)^f^'^^(x^) 



< Xf^'^^(x2) + (1-X)f^'^^(x^) 



(E.22) 



or 



CXf^(x2)+(l-X)f^(x^)][Xf(x2)+(l-X)f(x^)] < Xf^'*'l(x2)+(1-X)f^'^^(x^) . 

(E.23) 
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This Is the inequality we set out to show in (E.16). The 
theorem is proved. 

It has been shown that if f(x) is convex and f £ 0 then 
f (x) is convex. By Theorem 3> f (x) is convex. By 
Theorem 3 again, f (x) is convex. This reasoning can be 
followed for all powers K where K is a positive integer. 

The basic theorem is established. 



207 



LIST OF REFERENCES 



1. Bryson, A. E., "Applications of Optimal Control Theory 
in Aerospace Engineering," Journal of Spacecraft and 
Rockets , V. 4, no. 5, p. 5^5-553, May 1967. 

2. Bryson, A. E. and Denham, V/. F., "A Steepest-Ascent 
Method for Solving Optimum Programming Problems," 

Journal of Applied Mechanics , p. 2^7-257, June 1962. 

3. Bryson, A. E. and Hedrick, J. K., "Minimum Time Turns 
for a Supersonic Airplane at Constant Altitude," 

Journal of Aircraft , v. 8, no. 3, p. 182-187, March 1971. 

4. Bryson, A. E., Hoffman, W. C. , and Desai, M. N., The 
Energy-State Approximation in Performance Optimization 
of Supersonic Aircraft , AIAA Paper 68-877 » Pasadena , 
California, August 1968 . 

5. Kelley, H. J. and Edelbaum, T. N. , "Energy Climbs, Energy 
Turns, and Asymptotic Expansions," Journal of Aircraft , 

V. 7> no. 1, p. 93-95, January-February 1970. 

6. Kelley, H. J. and Lefton, L. , Differential Turns , paper 
presented at AIAA Atmospheric Flight Mechanics Special- 
ists Conference, Palo Alto, California, 11-13 September 
1972. 

7. Landgraf, S. K., "Some Applications of Performance 
Optimization Techniques to Aircraft," Journal of Aircraft, 
V. 2, no. 2, p. 153-154, March-April ISW. 

8. Taylor, L. W., Smith, H. J. , and Iliff, K. W., "Experience 
Using Balakrishnan ' s Epsilon Technique to Compute Optimum 
Flight Profiles," Journal of Aircraft , v. 7, no. 2, 

p. 182-197, March-April 1970. 

9. Fiacco, A. V. and McCormick, G. P., Nonlinear Programming: 
Sequential Unconstrained Minimization Techniques, V/iley, 

I96F: 

10. Jones, A. P. and McCormick, G. P., "A Generalization of 

the Method of Balakrishnan: Inequality Constraints and 

Initial Conditions," SIAM Journal on Control, p. 218-225, 
May 1970. 

11. Taylor, J. M. , Optimization: Application of the Epsilon 

Method , Ph.D. Thesis, Dept, of Electrical Engineering, 
University of Wyoming, Laramie, 1970. 



\ 



208 



12. Taylor, J. M. and Constantinides , C. T., "Optimization: 
Application of the Epsilon Method," IEEE Transactions on 
Automatic Control , p. 128-131, February 1972. 

13 . Department of Engineering, University of California 
Los Angeles Report 67-61, On a New Computing Technique 

in Optimal Control , by A. V. Balakrishnan, December I 967 . 

IM. Balakrishnan, A. V., On a Nevf Computing Technique in 
Optimal Control and Its Application to Minimal Time 
Flight Profile Optimization , paper presented at Inter- 
national Colloquium on Methods of Optimization, 2nd, 
Novosibirsk, USSR, June 1968. 

15 . Athans, M. and Falb, P. L. , Optimal Control: An 

Introduction to the Theory and Its Applications, 
McGraw-Hill, 1966. 

16 . Kirk, D. E., Optimal Control Theory: an Introduction , 
Prentice-Hall, 1970. 

17 . Rockafellar, R. T., Convex Analysis , Princeton, 1970. 

18 . Rockafellar, R. T., "Integrals Ifliich are Convex 
Functionals," Pacific Journal of Mathematics, v. 24, 
no. 3, P. 525 - 539 “, i968'. 

19 . Mlele, A., Flight Mechanics: Theory of Flight Paths , 

Addlson-Wesley , 1962. 



\ 



209 



t 



INITIAL DISTRIBUTION LIST 



No. 



1. Defense Documentation Center 
Cameron Station 
Alexandria! Virginia 2231^ 

2. Library, Code 0212 
Naval Postgraduate School 
Monterey, California 939^0 

3. Chairman, Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 939^0 

^1. Chairman, Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 939^0 

5. Professor D. E. Kirk, Code 52Ki 
Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 939^0 

6. Professor F. D. Faulkner, Code 53Fa 
Department of Mathematics 

Naval Postgraduate School 
Monterey, California 939^0 

7. Professor H. A. Titus, Code 52Ts 
Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 939^0 

8. Professor R. A. Hess, Code 57He 
Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 939^0 

9. Professor A.E. Fuhs, Code 57Fu 
Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 939^0 

10. Professor A. Gerba, Code 52Gz 
Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 939^0 



Copies 

2 

2 

1 

1 

3 

1 

1 

1 

1 

1 



210 



No. 



11. Dr. D. E. Zilmer, Code 607 
Naval Weapons Center 

China Lake, California 93555 

12. CDR. W. J. H. Smlthey 
Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 939^0 

13. CDR. Marie D. Hewett 
Naval Air Test Center 
Patuxent River, Maryland 



/ 



\ 



Copies 

1 

1 

3 



211 



SECURITY CLASSIFICATION OF THIS PAGE (Whmn Dmtm Entermd) 



REPORT DOCUMENTATIOH PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1. REPORT NUMBER 


2. GOVT ACCESSION NO. 


3. RECIPIENT’S CATALOG NUMBER 


4. title c«nd 5u6(ar.; 

A Second-Order Epsilon Method for 
Constrained Trajectory Optimization 


5. TYPE OF REPORT ft PERIOD COVERED 

Ph.D. Thesis; 

June, 197^ 


6. performing org. report number 


7. AuTMORf#; 

Marie David Hewett 


6. contract or grant number^*; 


9. PERFORMING organization NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 939^0 


10. program ELEMENT. PROJECT, TASK 
AREA 4 WORK UNIT NUMBERS 


n. CONTROLLING OFFICE NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 939^0 


12. REPORT DATE 

June 197^ 


13. NUMBER OF PAGES 

21i _ ^ 


14. monitoring agency NAME & ADDRESSC// different from Coniroittng Oiitce) 

Naval Postgraduate School 
Monterey, California 939^0 


15. SECURITY CLASS, (oi thte reporij 

Unclassified 


15«. DECLASSI FI CATION/ DOWN GRADING 
SCHEDULE 



16. DISTRIBUTION STATEMENT (o( thi9 Report) 

Approved for public release; distribution unlimited 



17. DISTRIBUTION STATEMENT (oi the ebetrmct entered in Biock 20, // different from Report) 



18. supplementary notes 



19. KEY WORDS (Continue on reveree eide ff neceeeejy end identify by biock number) 

optimization aircraft state 

trajectory limitation control 

Inequality penalty missile 

performance convexity 



20. ABSTRACT (Continue on reveree eide ff neceeeery end identify by biock number) 

A second-order epsilon method is developed for trajectory 
optimization problems. The method is applied to several 
aircraft and missile performance and air combat maneuvering 
problems. Heavy emphasis is placed on the realistic modeling 
of the flight vehicle’s motion and maneuvering limitations. 

The proposed optimization technique, which is an extension 
of Balakrlshnan ' s epsilon method, uses either the full 



DD 1473 



(Pago 1) 



EDITION OF 1 NOV 65 IS OBSOLETE 
S/N 0102-014- 6601 | 



SECURITY CLASSIFICATION OF THIS PAGE (f»ien Dete Sntered) 







I ^ ^ 



A 









1 







CfcCUKfTY CL ASSfFIC ATtON OF THIS PAGErw>i#n Dmtm Enfrmd) 



(20. continued) 

second-order Newton-Raphson method or the ’’modified" Hev;ton- 
Raphson method to minimize the epsilon functional. The full 
I'lewton-Raphson method exhibits terminal convergence characteristics 
superior to the ’’modified" method, vjhereas the ’’modified" method 
is generally superior in the Initial stages of a problem. An 
algorithm is developed which uses both techniques in a 
complementary way. 

A new penalty functional which has desirable theoretical 
properties and exhibits excellent computational behavior is 
introduced to treat state and control inequality constraints. 



DD Form 1473 
, 1 Jan 73 

S/N 0102-014-G601 



(BACK) 



213 



SECURITY CLASSIFICATION OF THIS P ACEfH7i*n Enfrmd) 





Thes i s 

H52624 

c.l 



Thests ■^•^^205 

H52624 Hewett 

A second-order epsi- 
lon method for con- 
strained trajectory 
optimization. 



' ri 9nc; 

U N> 



Hewett 

A second-order epsi- 
lon method for con- 
strained trajectory 
optimization. 



lhesH52624 



A second-order epsilon method for constr 




3 2768 001 91933 5 * 

DUDLEY KNOX LIBRARY C - / 



